4 Métodos lineales e ingenería de entradas

Una de las maneras más simples que podemos intentar para predecir \(y\) en función de las \(x_j\)´s es mediante una suma ponderada de los valores de las \(x_j's\), usando una función de la forma

\[f_\beta (x) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p,\] En este caso, escoger una función particular de esta familia, dada una muestra de entrenamiento \({\mathcal L}\), consiste en encontrar encontrar valores apropiados de las \(\beta\)’s, para construir un predictor:

\[\hat{f}(x) = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 \cdots + \hat{\beta} x_p\] y usaremos esta función \(\hat{f}\) para hacer predicciones \(\hat{y} =\hat{f}(x)\).

Ejemplos

Queremos predecir las ventas futuras anuales \(y\) de una tienda que se va a construir en un lugar dado. Las variables que describen el lugar son \(x_1 = trafico\_peatones\), \(x_2=trafico\_coches\). En una aproximación simple, podemos suponer que la tienda va a capturar una fracción de esos tráficos que se van a convertir en ventas. Quisieramos predecir con una función de la forma \[f_\beta (peatones, coches) = \beta_0 + \beta_1\, peatones + \beta_2\, coches.\] Por ejemplo, después de un análisis estimamos que

\(\hat{\beta}_0 = 1000000\) (ventas base, si observamos tráficos igual a 0: es lo que va a atraer la tienda)
\(\hat{\beta}_1 = (200)*0.02 = 4\) (capturamos 2% del tráfico peatonal, y cada capturado gasta 200 pesos)
\(\hat{\beta}_2 = (300)*0.01 =3\) (capturamos 1% del tráfico de autos, y cada capturado gasta 300 pesos)

Entonces haríamos predicciones con \[\hat{f}(peatones, coches) = 1000000 + 4\,peatones + 3\, coches.\] El modelo lineal es más flexible de lo que parece en una primera aproximación, porque tenemos libertad para construir las variables de entrada a partir de nuestros datos. Por ejemplo, si tenemos una tercera variable \(estacionamiento\) que vale 1 si hay un estacionamiento cerca o 0 si no lo hay, podríamos definir las variables

\(x_1= peatones\)
\(x_2 = coches\)
\(x_3 = estacionamiento\)
\(x_4 = coches*estacionamiento\)

Donde la idea de agregar \(x_4\) es que si hay estacionamiento entonces vamos a capturar una fracción adicional del trafico de coches, y la idea de \(x_3\) es que la tienda atraerá más nuevas visitas si hay un estacionamiento cerca. Buscamos ahora modelos de la forma \[f_\beta(x_1,x_2,x_3,x_4) = \beta_0 + \beta_1x_1 + \beta_2 x_2 + \beta_3 x_3 +\beta_4 x_4\] y podríamos obtener después de nuestra análisis las estimaciones - \(\hat{\beta}_0 = 800000\) (ventas base) - \(\hat{\beta}_1 = 4\) - \(\hat{\beta}_2 = (300)*0.005 = 1.5\) - \(\hat{\beta}_3 = 400000\) (ingreso adicional si hay estacionamiento por nuevo tráfico) - \(\hat{\beta}_4 = (300)*0.02 = 6\) (ingreso adicional por tráfico de coches si hay estacionamiento)

y entonces haríamos predicciones con el modelo \[\hat{f} (x_1,x_2,x_3,x_4) = 800000 + 4\, x_1 + 1.5 \,x_2 + 400000\, x_3 +6\, x_4\]

Variables derivadas

Nótese que cuando usamos modelos lineales, es necesario considerar qué entradas usar para construir el modelo de forma que la estructura lineal sea razonable. - Estas entradas son construidas a partir de las variables originales - Agregar variables tiene dos efectos: disminuye el sesgo, y también puede reducir el sesgo irreducible. - Omitir o simplificar estos modelos disminuyen típicamente la varianza pero pueden aumentar el sesgo irreducible.

En este caso particular, teníamos peatones, coches y estacionamiento, pero construimos también el producto de tráfico de coches y estacionamiento, pues la relación de ventas con coches es diferente dependiendo de sí la tienda tiene estacionamiento o no.

4.1 Aprendizaje de coeficientes (ajuste)

En el ejemplo anterior, los coeficientes fueron calculados (o estimados) usando experiencia, reglas, argumentos teóricos, o quizá otras fuentes de datos (como estudios o encuestas, conteos, etc.)

Ahora quisiéramos construir un algoritmo para aprender estos coeficientes del modelo \[f_\beta (x) = \beta_0 + \beta_1 x_1 + \cdots \beta_p x_p\] a partir de una muestra de entrenamiento de datos históricos de tiendas que hemos abierto antes: \[{\mathcal L}=\{ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}), \ldots, (x^{(N)}, y^{(N)}) \}\] El criterio de ajuste (algoritmo de aprendizaje) más usual para regresión lineal es el de mínimos cuadrados.

Construimos las predicciones (ajustados) para la muestra de entrenamiento: \[ f_\beta (x^{(i)}) = \beta_0 + \beta_1 x_1^{(i)}+ \cdots + \beta_p x_p^{(i)}\]

Y consideramos las diferencias de los ajustados con los valores observados:

\[e^{(i)} = y^{(i)} - f_\beta (x^{(i)})\]

La idea entonces es minimizar la suma de los residuales al cuadrado, para intentar que la función ajustada pase lo más cercana a los puntos de entrenamiento que sea posible. La función de pérdida que utilizamos más frecuentemente es la pérdida cuadrática, dada por:

\[L(\beta) = \frac{1}{N}\sum_{i=1}^N (y^{(i)} - f_\beta(x^{(i)}))^2\]

Mínimos cuadrados para regresión lineal

Buscamos encontrar: \[\hat{\beta} = \mathrm{arg\,min}_{\beta} L(\beta) = \mathrm{arg\,min}_{\beta}\frac{1}{N}\sum_{i=1}^N (y^{(i)} - f_\beta(x^{(i)}))^2\] donde \[f_\beta (x^{(i)}) = \beta_0 + \beta_1 x_1^{(i)}+ \cdots + \beta_p x_p^{(i)}\]

Hay varias maneras de resolver este problema: puede hacerse analíticamente con álgebra lineal, o con algún método numérico como descenso máximo (que puede escalarse fácilmente). Típicamente la función objetivo es convexa, y la solución es única, excepto en casos degenerados que podremos evitar más adelante usando regularización.

Observación: Como discutimos al final de la sección anterior, minimizar directamente el error de entrenamiento para encontrar los coeficientes puede resultar en en un modelo sobreajustado/con varianza alta/ruidoso. Hay cuatro grandes estrategias para mitigar este problema: restringir o estructurar la familia de funciones, penalizar la función objetivo, perturbar la muestra de entrenamiento, o cambiar el proceso de minimización perturbando la función objetivo en cada paso o deteniendo el proceso antes de llegar a un mínimo sobreajustado. El método mas común es cambiar la función objetivo, que discutiremos más adelante en la sección de regularización.

4.2 Ingeniería de entradas

Algunas veces, encontrar la estructura apropiada puede requerir más trabajo que simplemente escoger una familia de modelos. Por ejemplo, en el caso de precios de casa, vimos que podríamos mejorar el ajuste haciendo que el coeficiente de área habitable dependiera de la calidad de los terminados @ref(medicioncostosa).

Usualmente tendremos que hacer varias transformaciones para obtener buen desempeño de un modelo lineal. En la siguientes secciones mostramos algunas de las más usuales.

4.3 Variables categóricas

En primer lugar, podemos incluir variables categóricas creando variables numéricas 0-1 para cada categoría. Por ejemplo para la variable calidad sótano:

library(tidyverse)
library(tidymodels)
library(gt)
source("../R/casas_traducir_geo.R")
set.seed(83)
casas_split <- initial_split(casas, prop = 0.75)
casas_entrena <- training(casas_split)
casas_entrena |> count(calidad_sotano)

# A tibble: 5 × 2
  calidad_sotano     n
  <chr>          <int>
1 Ex                89
2 Fa                23
3 Gd               457
4 TA               500
5 <NA>              26

El mejor nivel es Ex (excelente), luego sigue Gd (bueno), luego Fa (razonable) y finalmente TA (típico)). Hay otro nivel Po (Malo) que no aparece en estos datos.

En primer lugar, podemos codificar los valores faltantes, que en este caso indican casas sin sótano:

receta_na <- recipe(~ calidad_sotano, casas_entrena) |> 
  step_unknown(calidad_sotano, new_level = "no_sótano") |> 
  step_relevel(calidad_sotano, ref_level = "TA")
casas_preproc <- prep(receta_na) |> juice()
casas_preproc |> count(calidad_sotano)

# A tibble: 5 × 2
  calidad_sotano     n
  <fct>          <int>
1 TA               500
2 Ex                89
3 Fa                23
4 Gd               457
5 no_sótano         26

Ahora convertimos a codificación dummy:

set.seed(7)
receta_dummy <- 
  recipe( ~ calidad_sotano, casas_entrena) |> 
  step_unknown(calidad_sotano, new_level = "no_sótano") |> 
  step_relevel(calidad_sotano, ref_level = "TA") |> 
  step_dummy(calidad_sotano, keep_original_cols = TRUE)
# preparar receta
receta_dummy_prep <- prep(receta_dummy) 
# extrae los datos de entrenamiento preprocesados
receta_dummy_prep |> juice() |> 
  sample_n(10) |> gt() |> 
  tab_options(table.font.size = 10)

calidad_sotano	calidad_sotano_Ex	calidad_sotano_Gd
TA	0	0
Gd	0	1
Ex	1	0
Ex	1	0
TA	0	0
TA	0	0
TA	0	0
TA	0	0
TA	0	0
TA	0	0

Nótese que no hay columna para el nivel TA, que tomamos como referencia. Incluir esta columna sería redundante, pues tenemos una constante en el predictor. En general, cuando una variable categórica tiene \(k\) niveles, esta codificación produce \(k-1\) columnas binarias.

Veamos qué pasa cuando preprocesamos datos de prueba (para después poder hacer predicciones):

prueba_casas <- testing(casas_split)
# supongamos que un nuevo nivel aparece
prueba_casas$calidad_sotano[1] <- "no visto antes"
datos <- bake(receta_dummy_prep, prueba_casas)

Warning: There are new levels in a factor: no visto antes
New levels will be coerced to `NA` by `step_unknown()`.
Consider using `step_novel()` before `step_unknown()`.

Warning: There are new levels in a factor: NA

En este caso, podemos hacer nuestro flujo más robusto incluyendo un nuevo nivel en los factores donde pondremos casos no vistos. Modificamos nuestra receta:

receta_dummy <- 
  recipe( ~ calidad_sotano, casas_entrena) |> 
  step_novel(calidad_sotano, new_level = "nuevo") |> 
  step_unknown(calidad_sotano, new_level = "no_sótano") |> 
  step_relevel(calidad_sotano, ref_level = "TA") |> 
  step_dummy(calidad_sotano, keep_original_cols = TRUE)
# preparar receta
receta_dummy_prep <- prep(receta_dummy)

prueba_casas <- testing(casas_split)
# supongamos que un nuevo nivel aparece
prueba_casas$calidad_sotano[1] <- "no visto antes"
datos <- bake(receta_dummy_prep, prueba_casas)
datos |> head() |> gt() |> tab_options(table.font.size = 10)

calidad_sotano	calidad_sotano_Ex	calidad_sotano_Gd	calidad_sotano_nuevo
nuevo	0	0	1
Ex	1	0	0
TA	0	0	0
Gd	0	1	0
TA	0	0	0
TA	0	0	0

Y podemos ignorar el nuevo nivel al hacer predicciones (que equivale a ponerlo en la categoría de referencia, que en este caso es TA), o podemos lidiar de manera ad-hoc con este nivel.

Otro problema con el que podemos encontrarnos es variables categóricas que son muy ralas. Por ejemplo, una variable que tiene muchas categorías y algunas de ellas tienen muy pocos datos, además de que es probable que observemos nuevas categorías en el futuro. Por ejemplo, para la variable de zona:

casas_entrena |> count(nombre_zona) |> 
  arrange(desc(n)) |> gt() |> tab_options(table.font.size = 10)

nombre_zona	n
NAmes	163
CollgCr	113
OldTown	80
Edwards	74
Somerst	65
Gilbert	61
Sawyer	59
NridgHt	58
NWAmes	57
BrkSide	45
SawyerW	42
Crawfor	39
Mitchel	36
IDOTRR	31
NoRidge	30
Timber	26
ClearCr	22
StoneBr	20
SWISU	18
Blmngtn	14
BrDale	14
MeadowV	12
NPkVill	8
Veenker	7
Blueste	1

En este caso, tenemos muchas categorías, algunas con muy pocos datos, y es posible que observemos nuevos datos. Una técnica es agrupar los datos de baja cardinalidad en un nuevo nivel (incluyendo categorías no observadas en entrenamiento):

receta_vecindario_1 <- 
  recipe( ~ nombre_zona, casas_entrena) |>
  step_other(nombre_zona, threshold = 0.01, other = "otras") 
receta_vecindario <- receta_vecindario_1 |> 
  step_dummy(nombre_zona)
# preparar receta
receta_vecindario_prep <- prep(receta_vecindario_1)
set.seed(8231)
receta_vecindario_prep |> juice() |> 
  count(nombre_zona) |> arrange(desc(n)) |> gt()

nombre_zona	n
NAmes	163
CollgCr	113
OldTown	80
Edwards	74
Somerst	65
Gilbert	61
Sawyer	59
NridgHt	58
NWAmes	57
BrkSide	45
SawyerW	42
Crawfor	39
Mitchel	36
IDOTRR	31
NoRidge	30
Timber	26
ClearCr	22
StoneBr	20
SWISU	18
otras	16
Blmngtn	14
BrDale	14
MeadowV	12

En este caso, las zonas de baja frecuencia fueron agrupadas en la categoría “otras”. Si observamos un nuevo nivel al momento de predicción:

prueba_casas <- testing(casas_split)
# supongamos que un nuevo nivel aparece
prueba_casas$nombre_zona[1] <- "Xochimilco"
datos <- bake(prep(receta_vecindario), prueba_casas)
datos |> head() |>  gt() |> tab_options(table.font.size = 10)

nombre_zona_CollgCr	nombre_zona_NAmes	nombre_zona_Sawyer	nombre_zona_Somerst	nombre_zona_otras
0	0	0	0	1
0	0	0	1	0
0	0	1	0	0
1	0	0	0	0
0	1	0	0	0
0	1	0	0	0

El proceso general es (ver por ejemplo esta lista):

Variables categóricas

Establecemos los niveles que puede tener cada variable, incluyendo la posibilidad de categorías nuevas al momento de predecir, y categorías para valores no disponibles (NAs) (es posible también imputar con algún método en caso necesario).
Reorganizamos factores dependiendo del problema. Por ejemplo, incluir categorías de baja frecuencia en una categoría separada, o manipulaciones ad-hoc dependiendo del problema.
Sustituimos variables categóricas con \(K\) niveles en \(K-1\) columnas indicadoras de los niveles (estableciendo) alguna categoría como referencia. Esto no es estrictamente necesario en otros métodos, o si utilizamos regularización (ver sección siguiente).

4.4 Interacciones

Otra manera de expandir nuestro modelo es la utilización de interacciones, que muchas veces son clave para tener éxito con modelos lineales. Vimos ejemplos de interacciones en el ejemplo de las casas (@ref(medicioncostosa)) y en el primer ejemplo de ventas de tiendas que dependían del tráfico.

Ejemplo

Si \(x_1\) es el área en metros cuadrados de una casa, y \(x_2\) una calificación numérica de su calidad, podemos considerar el modelo sin interacciones:

\[\beta_0 + \beta_1x_1 + \beta_2 x_2\]

Pero no tiene mucho sentido que el efecto marginal de \(x_1\) sea constante para cualquier nivel de calidad, y tampoco que la calidad de terminados agregue una cantidad fija al precio de la casa sin tomar en cuenta su tamaño. Podemos remediar esto creando una nueva variable que es le producto de \(x_1\) y \(x_2\):

\[x_3 = x_1 x_2\]

y agregando, nuestro predictor para precio es

\[ \beta_0 + \beta_1x_1 + \beta_2 x_2 + \beta_3 x_1x_2\]

Ahora notemos que para \(x_2\) fija, el modelo es

\[(\beta_0 + \beta_2x_2) + (\beta_1 + \beta_3x_2)x_1 = \gamma_0 + \gamma_1x_1\] De modo que es lineal en \(x_1\). La diferencia es que cuando cambia \(x_2\), la recta que ajustamos es diferente.

beta <- c(0, 50, 100, 20)
combs_tbl <- crossing(x_1 = seq(2, 20, by = 1), x_2 = seq(0, 10, by = 2)) |> 
  mutate(x_3 = x_1 * x_2) |> 
  mutate(pred = beta[1] + beta[2]*x_1 + beta[3]*x_2 + beta[4]*x_3)
ggplot(combs_tbl, aes(x = x_1, y = pred, group = x_2, colour = x_2)) +
  geom_line()

Y vemos que cuando la calidad es baja, el precio por metro cuadrado es más bajo que cuando la calidad es alta. Otra manera de pensar esto es que la inclusión de la interacción produce curvas marginales que rotan dependiendo del valor de otras variables.

Pregunta. ¿puedes pensar en otros casos donde las interacciones deben jugar un papel importante?

Interacciones

Transformamos las variables categóricas a dummies. Transformamos las variables numéricas si es necesario (normalizar, aplicar transformación no lineal, etc.)
Incluimos interacciones de la siguiente forma:

Para la interaccion de dos variables numéricas \(x_1\) y \(x_2\) agregamos el producto \(x_3 = x_1x_2\).
Para interacción de una variable categórica \(g\) con una numérica \(x\) podemos hacer el mismo procedimiento multiplicando la variable categórica por cada una de las variable dummy que creamos a partir de \(g\). Esto en efecto produce una pendiente para \(x\) dependiendo del valor que toma \(g\).

4.5 Ejemplo: precios de casas

En el ejemplo de precios de casas, por ejemplo, es claro que el efecto en ventas del tamaño de las áreas (habitable, garage, etc.) depende de la calidad de los terminados, como vimos en la introducción. En la siguiente receta de preprocesamiento:

Cortamos calidad general en 5 grupos: este paso no es necesario y puede dañar el desempeño, pero es consistente con el análisis que hicimos anteriormente.
Lidiamos con niveles nuevos y los ponemos en una categoría “nuevo” (para que nuestro modelo no falle al momento de predicción)
Ponemos los faltantes de calidad sotano y garage en una categoría nueva (no tienen sótano y/o garage)
Agrupamos las zonas con pocas observaciones en una categoría de “Otros”
Quitamos los NA’s de área garage y área sotano, que deben ser igual a 0 cuando no existen estas características.
Creamos variables dummy de todas las variables categóricas
Incluimos interacciones de distintas áreas con las dummy correspondientes, incluyendo zona con área habitable
Finalmente, eliminamos para el ajuste aquellas variables que tengan varianza cercana a cero (500 /1 quiere decir que elimina cualquier variable cuyo conteo del valor más común entre el conteo de la siguiente es mayor a 500).

library(tidyverse)
library(tidymodels)
library(gt)
source("../R/casas_traducir_geo.R")
set.seed(83)
casas_split <- initial_split(casas, prop = 0.75)
casas_entrena <- training(casas_split)
receta_casas <- recipe(precio_miles ~ 
           nombre_zona + 
           area_hab_m2 + area_garage_m2 + area_sotano_m2 + 
           area_lote_m2 + 
           año_construccion + 
           calidad_gral + calidad_garage + calidad_sotano + 
           num_coches  + 
           aire_acondicionado + condicion_venta, 
           data = casas_entrena) |> 
  step_filter(condicion_venta == "Normal") |> 
  step_select(-condicion_venta, skip = TRUE) |> 
  step_cut(calidad_gral, breaks = c(3, 5, 7, 8), 
           include_outside_range = TRUE) |>
  step_novel(nombre_zona, calidad_sotano, calidad_garage) |> 
  step_unknown(calidad_sotano, calidad_garage) |> 
  step_other(nombre_zona, threshold = 0.02, other = "otras") |> 
  step_mutate(area_sotano_m2 = ifelse(is.na(area_sotano_m2), 0, area_sotano_m2)) |> 
  step_mutate(area_garage_m2 = ifelse(is.na(area_garage_m2), 0, area_garage_m2)) |> 
  step_dummy(nombre_zona, calidad_gral, calidad_garage, calidad_sotano, aire_acondicionado) |> 
  step_interact(terms = ~ area_hab_m2:starts_with("calidad_gral")) |> 
  step_interact(terms = ~ area_hab_m2:starts_with("nombre_zona")) |> 
  step_interact(terms = ~ area_garage_m2:starts_with("calidad_garage")) |> 
  step_interact(terms = ~ area_sotano_m2: starts_with("calidad_sotano")) |> 
  step_nzv(all_predictors(), freq_cut = 500 / 1, unique_cut = 1)

Entrenamos la receta y vemos cuántos casos y columnas tenemos:

receta_casas_prep <- prep(receta_casas, verbose = TRUE)

oper 1 step filter [training] 
oper 2 step select [training] 
oper 3 step cut [training] 
oper 4 step novel [training] 
oper 5 step unknown [training] 
oper 6 step other [training] 
oper 7 step mutate [training] 
oper 8 step mutate [training] 
oper 9 step dummy [training] 
oper 10 step interact [training] 
oper 11 step interact [training] 
oper 12 step interact [training] 
oper 13 step interact [training] 
oper 14 step nzv [training] 
The retained training set is ~ 0.46 Mb  in memory.

datos_tbl <- juice(receta_casas_prep)
dim(datos_tbl)

[1] 907  62

datos_tbl |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  head() |> 
  gt()

area_hab_m2	area_garage_m2	area_sotano_m2	area_lote_m2	año_construccion	num_coches	precio_miles	nombre_zona_CollgCr	nombre_zona_Gilbert	nombre_zona_NAmes	nombre_zona_otras	calidad_gral_X.3.5.	calidad_gral_X.5.7.	calidad_gral_X.7.8.	calidad_garage_TA	calidad_sotano_Gd	calidad_sotano_TA	calidad_sotano_unknown	aire_acondicionado_Y	area_hab_m2_x_calidad_gral_X.3.5.	area_hab_m2_x_calidad_gral_X.5.7.	area_hab_m2_x_calidad_gral_X.7.8.	area_hab_m2_x_nombre_zona_CollgCr	area_hab_m2_x_nombre_zona_Gilbert	area_hab_m2_x_nombre_zona_NAmes	area_hab_m2_x_nombre_zona_otras	area_garage_m2_x_calidad_garage_TA	area_sotano_m2_x_calidad_sotano_Gd	area_sotano_m2_x_calidad_sotano_TA
137.22	20.07	79.99	584.55	1928	1	145.0	0	0	0	1	0	1	0	1	0	1	0	1	0.00	137.22	0.00	0.00	0.00	0.00	137.22	20.07	0.00	79.99
179.86	46.82	130.62	1039.96	2000	2	230.5	0	1	0	0	0	0	1	1	1	0	0	1	0.00	0.00	179.86	0.00	179.86	0.00	0.00	46.82	130.62	0.00
81.94	27.31	0.00	774.72	1959	1	106.5	0	0	1	0	1	0	0	1	0	0	1	1	81.94	0.00	0.00	0.00	0.00	81.94	0.00	27.31	0.00	0.00
153.01	20.07	74.88	637.13	1915	1	128.0	0	0	0	1	0	1	0	1	1	0	0	0	0.00	153.01	0.00	0.00	0.00	0.00	153.01	20.07	74.88	0.00
101.45	26.57	50.73	205.97	1970	1	88.0	0	0	0	1	1	0	0	1	0	1	0	1	101.45	0.00	0.00	0.00	0.00	0.00	101.45	26.57	0.00	50.73
71.35	36.79	71.35	668.90	1972	1	133.9	1	0	0	0	1	0	0	1	1	0	0	1	71.35	0.00	0.00	71.35	0.00	0.00	0.00	36.79	71.35	0.00

Finalmente, usamos un modelo lineal con las 62 entradas que acabamos de crear:

flujo_casas <- workflow() |> 
  add_recipe(receta_casas) |> 
  add_model(linear_reg() |> set_engine("lm"))
ajuste <- fit(flujo_casas, casas_entrena)

Aunque no es de interés particular para nosotros por el momento, examinamos los coeficientes (que no son tan simples de interpretar como discutiremos más adelante):

ajuste |> broom::tidy() |> 
  mutate(across(where(is.numeric), \(x) round(x, 2))) |> 
  select(term, estimate) |> 
  gt()

term	estimate
(Intercept)	-858.68
area_hab_m2	0.80
area_garage_m2	1.02
area_sotano_m2	0.88
area_lote_m2	0.01
año_construccion	0.39
num_coches	2.11
nombre_zona_CollgCr	16.70
nombre_zona_Crawfor	37.11
nombre_zona_Edwards	28.82
nombre_zona_Gilbert	31.31
nombre_zona_IDOTRR	11.81
nombre_zona_Mitchel	32.01
nombre_zona_NAmes	28.74
nombre_zona_NoRidge	22.00
nombre_zona_NridgHt	-6.51
nombre_zona_NWAmes	17.67
nombre_zona_OldTown	25.28
nombre_zona_Sawyer	37.91
nombre_zona_SawyerW	-11.95
nombre_zona_Somerst	-8.99
nombre_zona_Timber	3.23
nombre_zona_otras	-14.27
calidad_gral_X.3.5.	35.54
calidad_gral_X.5.7.	13.76
calidad_gral_X.7.8.	16.99
calidad_gral_X.8.max.	-67.26
calidad_garage_Fa	15.54
calidad_garage_Gd	39.04
calidad_garage_TA	18.71
calidad_garage_unknown	20.61
calidad_sotano_Fa	72.92
calidad_sotano_Gd	54.95
calidad_sotano_TA	60.56
calidad_sotano_unknown	59.88
aire_acondicionado_Y	15.17
area_hab_m2_x_calidad_gral_X.3.5.	-0.25
area_hab_m2_x_calidad_gral_X.5.7.	0.05
area_hab_m2_x_calidad_gral_X.7.8.	0.19
area_hab_m2_x_calidad_gral_X.8.max.	0.81
area_hab_m2_x_nombre_zona_CollgCr	-0.25
area_hab_m2_x_nombre_zona_Crawfor	-0.14
area_hab_m2_x_nombre_zona_Edwards	-0.39
area_hab_m2_x_nombre_zona_Gilbert	-0.33
area_hab_m2_x_nombre_zona_IDOTRR	-0.21
area_hab_m2_x_nombre_zona_Mitchel	-0.42
area_hab_m2_x_nombre_zona_NAmes	-0.31
area_hab_m2_x_nombre_zona_NoRidge	-0.16
area_hab_m2_x_nombre_zona_NridgHt	-0.03
area_hab_m2_x_nombre_zona_NWAmes	-0.26
area_hab_m2_x_nombre_zona_OldTown	-0.32
area_hab_m2_x_nombre_zona_Sawyer	-0.43
area_hab_m2_x_nombre_zona_SawyerW	-0.07
area_hab_m2_x_nombre_zona_Somerst	0.00
area_hab_m2_x_nombre_zona_Timber	-0.14
area_hab_m2_x_nombre_zona_otras	0.00
area_garage_m2_x_calidad_garage_Fa	-0.59
area_garage_m2_x_calidad_garage_Gd	-0.97
area_garage_m2_x_calidad_garage_TA	-0.78
area_sotano_m2_x_calidad_sotano_Fa	-0.86
area_sotano_m2_x_calidad_sotano_Gd	-0.53
area_sotano_m2_x_calidad_sotano_TA	-0.67

Nótese que:

En esta tabla están los coeficientes \(\beta_i\) en las covariables que creamos a partir de las variables de entrada.
El modelo lineal no tiene que ser lineal en las variables que recibimos originalmente en la tabla de datos.
En este ejemplo, convertimos algunas variables a dummy, y multiplicamos algunas variables de área por esas variables dummy.

Finalmente, evaluamos el desempeño sobre las ventas normales:

metricas <- metric_set(mape, mae, rmse, rsq)
casas_prueba_normal <- testing(casas_split) |> 
  filter(condicion_venta == "Normal")
metricas(casas_prueba_normal |> bind_cols(predict(ajuste, casas_prueba_normal)), 
     truth = precio_miles, estimate = .pred) |> 
  gt() |> fmt_number(.estimate, decimals = 2)

.metric	.estimator	.estimate
mape	standard	10.34
mae	standard	17.12
rmse	standard	23.49
rsq	standard	0.89

4.6 No linealidad y atípicos

En un primer ejemplo consideremos la variable de área de lote:

library(patchwork)
g_1 <- ggplot(casas_entrena, aes(x = area_lote_m2, y = precio_miles)) +
  geom_point() + #geom_smooth(method = "loess", span = 0.5, se = FALSE,
                #             method.args = list(degree = 1)) +
  geom_smooth(method = "lm", se = FALSE)
g_2 <- ggplot(casas_entrena |> filter(area_lote_m2 < 5000), aes(x = area_lote_m2, y = precio_miles)) +
  geom_point(data = casas_entrena, aes(colour = area_lote_m2 > 5000)) + 
  #geom_smooth(method = "loess", span = 0.5, se = FALSE,
  #                           method.args = list(degree = 1)) +
  geom_smooth(method = "lm", se = FALSE)
g_1 + g_2

Y notamos que hay algunos valores grandes que pueden perturbar el ajuste lineal. Esto puede producir varianza alta en las predicciones, pues el ajuste depende mucho de unos cuantos valores de entrenamiento. Una solución puede ser transformar la entrada por ejempo usando el logaritmo, que comprime la cola derecha de la distribución de la variable que tiene mucho sesgo:

ggplot(casas_entrena, aes(x = area_lote_m2, y = precio_miles)) +
  geom_point() +
  geom_smooth(method = "loess", span = 0.5, se = FALSE) +
  scale_x_log10()

`geom_smooth()` using formula = 'y ~ x'

Nótese que probablemente tendremos que agregar más flexibilidad en nuestro predictor para capturar apropiadamente la información en esta variable.

4.7 No linealidad y splines

En algunos casos, la relación de una variable de entrada con la predicción es no lineal. Podemos entonces incluír entradas derivadas de la original usando transformaciones no lineales: por ejemplo, transformar entradas usando el logaritmo, o agregar el cuadrado o la raíz de las variables de entrada.

Una de las maneras más simples y menos problemáticas de hacer esto es usando splines naturales para modelar, que son funciones cúbicas por tramos dos veces diferenciables. Los tramos están definidos por nudos que podemos definir por ejemplo igualmente espaciados en los datos.

valores_x <- seq(-10, 110, 1)
base_splines <- splines::ns(x = valores_x, 
  knots = c(33, 66), Boundary.knots = c(0, 100))
spline_1 <- base_splines %*% c(1, 1, 1)
spline_2 <- base_splines %*% c(-0.1, 2, 1)
tibble(x = valores_x, y = spline_1, spline = 1) |> 
bind_rows(tibble(x = valores_x, y = spline_2, spline = 2)) |> 
  ggplot(aes(x = x, y = y, 
    group = spline, colour = factor(spline))) + 
  geom_point() + geom_line() +
  geom_vline(xintercept = c(0, 33, 66, 100), 
    colour = "red")

Estos dos son ejemplos de funciones cúbicas por tramos y dos veces diferenciables, con nudos en 0, 33, 66 y 100. Su forma particular depende de tres coeficientes, que pueden pensarse también como definidos por dónde tienen que pasar la curva en \(y\) para los valores \(x = 33, 66\) y \(100\). Extrapolan linealmente fuera del rango de los datos.

La ventaja de utilizar estos splines es que son estables en el cálculo, pues a lo más utilizan potencias cúbicas, y la complejidad puede aumentarse incrementando el número de nodos.

Ejemplo

Revisamos nuestro ejemplo de rendimiento de coches:

library(tidyverse)
library(gt)
auto <- read_csv("../datos/auto.csv")
datos <- auto[, c('name', 'weight','year', 'mpg', 'displacement')]
datos <- datos |> mutate(
  peso_kg = weight * 0.45359237,
  rendimiento_kpl = mpg * (1.609344 / 3.78541178), 
  año = year)

Vamos a separar en muestra de entrenamiento y de prueba estos datos. Podemos hacerlo como sigue (75% para entrenamiento aproximadamente en este caso, así obtenemos alrededor de 100 casos para prueba):

library(tidymodels)
set.seed(121)
datos_split <- initial_split(datos, prop = 0.75)
datos_entrena <- training(datos_split)
datos_prueba <- testing(datos_split)

Vamos a usar año y peso de los coches para predecir su rendimiento:

ggplot(datos_entrena, aes(x = peso_kg, y = rendimiento_kpl, colour = año)) +
  geom_point()

Nuestra receta incluye la transformación no lineal de splines:

receta_lineal <- recipe(rendimiento_kpl ~ peso_kg + año, datos_entrena) |> 
  step_ns(peso_kg, deg_free = 3) |> 
  step_ns(año, deg_free = 2)
mod_lineal <- linear_reg() |>  
  set_engine("lm")  
flujo <- workflow() |>  
  add_recipe(receta_lineal) |> 
  add_model(mod_lineal)

Los datos de entrada son los siguientes:

juice(prep(receta_lineal)) |> head() |> gt()

rendimiento_kpl	peso_kg_ns_1	peso_kg_ns_2	peso_kg_ns_3	año_ns_1	año_ns_2
12.754311	-0.1384501	0.3919536	-0.2338715	0.1263158	-0.08343893
13.136941	-0.1508427	0.4978146	-0.2970367	0.5652102	0.00590925
12.754311	0.3794499	0.4646633	-0.2232209	0.4765544	0.35513660
7.695101	0.4471996	0.4245669	-0.1625354	0.5652102	0.00590925
8.757960	0.4368600	0.4313780	-0.1743974	0.5652102	0.00590925
14.242314	-0.1126420	0.2985773	-0.1781555	0.5794245	-0.12316573

Nótese que tenemos 4 entradas en lugar de las 2 originales, pues creamos dos transformaciones no lineales de peso_kg. El modelo es lineal en estas 4 variables, pero no en las 2 originales. Ajustamos:

flujo_ajustado <- fit(flujo, datos_entrena)

Y ahora podemos graficar los resultados y vemos cómo pudimos capturar la relación no lineal entre peso y rendimiento:

dat_graf <- tibble(peso_kg = seq(900, 2200, by = 10)) |>   
  crossing(tibble(año = c(70, 75, 80)))
dat_graf <- dat_graf |> 
  mutate(pred_1 = predict(flujo_ajustado, dat_graf) |> pull(.pred))
ggplot(datos_entrena, aes(x = peso_kg, group = año, colour = año)) +
  geom_point(aes(y = rendimiento_kpl), alpha = 0.6) + 
  geom_line(data = dat_graf, aes(y = pred_1),  size = 1.2)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Los grados de libertad también pueden afinarse utilizando un conjunto de validación como hicimos antes en vecinos más cercanos.

Ejemplo: casas

Por ejemplo, podríamos incluir un efecto no lineal de area_lote, calidad y condición general y año de construcción:

receta_casas <- recipe(precio_miles ~ 
           nombre_zona + 
           area_hab_m2 + area_garage_m2 + area_sotano_m2 + 
           area_lote_m2 + 
           año_construccion + 
           calidad_gral + calidad_garage + calidad_sotano + 
           condicion_gral + 
           num_coches  + 
           aire_acondicionado + condicion_venta, 
           data = casas_entrena) |> 
  step_filter(condicion_venta == "Normal") |> 
  step_select(-condicion_venta, skip = TRUE) |> 
  
  step_novel(nombre_zona, calidad_sotano, calidad_garage) |> 
  step_unknown(calidad_sotano, calidad_garage) |> 
  step_other(nombre_zona, threshold = 0.02, other = "otras") |> 
  step_mutate(area_sotano_m2 = 
    ifelse(is.na(area_sotano_m2), 0, area_sotano_m2)) |> 
  step_mutate(area_garage_m2 = 
    ifelse(is.na(area_garage_m2), 0, area_garage_m2)) |>
 # step_log(area_lote_m2) |> 
  step_ns(año_construccion, deg_free = 2) |> 
  step_ns(calidad_gral, deg_free = 2) |> 
  step_ns(condicion_gral, deg_free = 2) |> 
  step_ns(area_lote_m2, deg_free = 3) |> 
  step_dummy(nombre_zona,  calidad_garage, calidad_sotano, aire_acondicionado) |> 
  step_interact(terms = ~ area_hab_m2:starts_with("calidad_gral")) |> 
  step_interact(terms = ~ area_hab_m2:starts_with("nombre_zona")) |> 
  step_interact(terms = ~ area_garage_m2:starts_with("calidad_garage")) |> 
  step_interact(terms = ~ area_sotano_m2: starts_with("calidad_sotano")) |> 
  step_nzv(all_predictors(), freq_cut = 500 / 1, unique_cut = 1)

flujo_casas <- workflow() |> 
  add_recipe(receta_casas) |> 
  add_model(linear_reg() |> set_engine("lm"))
ajuste <- fit(flujo_casas, casas_entrena)

Finalmente, evaluamos el desempeño sobre las ventas normales, y obtenemos una mejoría con respecto a nuestro modelo anterior:

metricas <- metric_set(mape, mae, rmse, rsq)
casas_prueba_normal <- testing(casas_split) |> 
  filter(condicion_venta == "Normal")
metricas(casas_prueba_normal |> 
  bind_cols(predict(ajuste, casas_prueba_normal)), 
     truth = precio_miles, estimate = .pred) |> 
  gt() |> fmt_number(.estimate, decimals = 2)

.metric	.estimator	.estimate
mape	standard	8.47
mae	standard	14.45
rmse	standard	19.77
rsq	standard	0.92

Finalmente, examinamos la respuesta contra la predicción:

casas_prueba_normal |> 
  bind_cols(predict(ajuste, casas_prueba_normal)) |> 
  ggplot(aes(x = .pred, y = precio_miles)) + geom_abline() + 
    geom_point(colour = "red") + coord_obs_pred()