paquetes <- c("readxl", "dplyr", "tidyr", "ggplot2", "knitr", "kableExtra")
faltan <- paquetes[!(paquetes %in% installed.packages()[, "Package"])]
if (length(faltan) > 0) install.packages(faltan)

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
library(kableExtra)

theme_set(theme_minimal(base_size = 12))

# para escribir los valores p en el texto
valor_p <- function(p) ifelse(p < 0.001, "p < 0.001", paste("p =", round(p, 3)))

Contexto y problema

Elegí la base 1, que contiene 314 anuncios de Mazda2 publicados en OLX Colombia. La pregunta que quiero responder es qué tan bien se puede estimar el precio anunciado de un Mazda2 usado a partir de sus características (año del modelo, kilometraje y tipo de transmisión, principalmente).

Este tipo de estimación le sirve a quien va a vender su carro y necesita fijar un precio de publicación, a un comprador que quiere saber si un anuncio está caro frente a carros parecidos, o a una compraventa que recibe un carro en parte de pago y quiere una primera referencia antes del peritaje.

Hay que tener presente desde el inicio que la variable es el precio que pide el vendedor, no el precio al que finalmente se vende el carro. Los anuncios no dicen nada sobre la negociación, así que todo lo que se concluya aplica al precio de publicación.

1.1 Descripción de las variables y calidad de los datos

Carga de la base

mazda <- read_excel("base1_mazdacol.xlsx")
mazda$id <- 1:nrow(mazda)
dim(mazda)
## [1] 314  12
mazda %>%
  select(precio, kilometraje, transmision, modelo, color, ciudad, departamento) %>%
  head() %>%
  kable(caption = "Primeras filas de la base") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Primeras filas de la base
precio kilometraje transmision modelo color ciudad departamento
34000000 140000 Automática 2008 Gris Medellín Antioquia
26400000 112000 Automática 2012 Gris Medellín Antioquia
42500000 102222 Automática 2013 Gris Itagüí Antioquia
91000000 1000 Automática 2022 Gris Envigado Antioquia
35000000 134500 Automática 2012 Blanco Medellín Antioquia
67000000 61997 Automática 2017 Gris Medellín Antioquia 
# fecha de extracción (la columna web-scraper-order empieza con un timestamp)
marca <- as.numeric(sub("-.*", "", mazda$`web-scraper-order`))
fecha_extraccion <- as.Date(as.POSIXct(min(marca), origin = "1970-01-01"))
fecha_extraccion
## [1] "2022-09-12"

La base tiene 314 filas. Las variables útiles para el análisis son:

  • precio: precio anunciado en pesos (variable respuesta). Lo trabajo en millones para que los coeficientes sean fáciles de leer.
  • kilometraje: kilómetros reportados por el vendedor.
  • transmision: Automática o Mecánica (categórica binaria).
  • modelo: año del modelo. Como los datos se extrajeron en septiembre de 2022 (2022-09-12), calculé la edad del carro como 2022 menos el año del modelo.
  • color, ciudad y departamento: categóricas con muchos niveles.
  • link: aunque se llama así, contiene el texto del anuncio (precio, año, km, título y barrio). No se usa como predictor, pero sirve para revisar la calidad de los datos.

Las columnas web-scraper-order, web-scraper-start-url y link-href son información del proceso de extracción.

Valores faltantes y duplicados

vars <- c("precio", "kilometraje", "transmision", "modelo", "color", "ciudad", "departamento")
colSums(is.na(mazda[, vars]))
##       precio  kilometraje  transmision       modelo        color       ciudad 
##            0            0            0            0            0            0 
## departamento 
##            0
# anuncios repetidos en todas las variables de análisis
sum(duplicated(mazda[, vars]))
## [1] 1
# enlaces repetidos
sum(duplicated(mazda$`link-href`))
## [1] 0

No hay valores faltantes y solo aparece 1 anuncio duplicado. Sin embargo, al leer el texto de los anuncios encontré problemas que no se ven en las columnas numéricas:

  1. Hay anuncios que no son de un Mazda2: CX-5, CX-9, CX-30, Mazda 3, Mazda 6, una camioneta de estacas de 1995, una “camioneta Mazda Touring” y dos Renault Stepway que ponen “Mazda 2” en el título solo para aparecer en la búsqueda.
  2. Un mismo Mazda2 2008 de salvamento (declarado pérdida total por una aseguradora) está publicado tres veces. Un carro de salvamento tiene un mercado distinto, así que no es comparable con el resto.
  3. Un Mazda2 2022 aparece a 6,8 millones, lo cual no es posible (probablemente le falta un dígito o es un anuncio engañoso).
  4. Un carro 2009 tiene 0 km, es decir, no se reportó el kilometraje.
  5. Once carros usados (modelo 2020 o anterior) tienen entre 35 y 183 km. Por el año y el precio es claro que el vendedor escribió el kilometraje en miles (por ejemplo, “116” en vez de 116.000).
  6. El anuncio con id 299 dice modelo 2019 en el campo del año, pero el título dice “MAZDA 2 2009”, y su precio (25,8 millones) y su kilometraje (129.000) corresponden a un 2009.

Los casos 1 a 4 se eliminan. Los casos 5 y 6 se corrigen, porque la información correcta se puede recuperar del mismo anuncio.

mazda$duplicado <- duplicated(mazda[, vars])

mazda <- mazda %>%
  mutate(titulo = tolower(link),
         otro_modelo = grepl("cx-?[0-9]|mazda ?3|mazda ?6|estacas|stepway|camioneta mazda", titulo),
         salvamento = grepl("salvamen", titulo),
         km_en_miles = modelo <= 2020 & kilometraje > 0 & kilometraje < 1000,
         motivo = case_when(
           duplicado ~ "Duplicado",
           otro_modelo ~ "No es un Mazda2",
           salvamento ~ "Salvamento",
           precio < 15000000 ~ "Precio imposible",
           modelo <= 2020 & kilometraje == 0 ~ "Km = 0 en carro usado",
           TRUE ~ "Se conserva"
         ),
         anuncio = substr(sub("^[^k]*km", "", link), 1, 60))

table(mazda$motivo)
## 
##             Duplicado Km = 0 en carro usado       No es un Mazda2 
##                     1                     1                    10 
##      Precio imposible            Salvamento           Se conserva 
##                     1                     2                   299
mazda %>%
  filter(motivo != "Se conserva") %>%
  transmute(id, precio_millones = precio / 1e6, modelo, kilometraje, motivo, anuncio) %>%
  kable(caption = "Anuncios eliminados") %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) %>%
  scroll_box(height = "300px")
Anuncios eliminados
id precio_millones modelo kilometraje motivo anuncio
2 26.4 2012 112000 No es un Mazda2 Vencambio Renault STEPWAY automático RECIBO TWINGO KIA PICAN
9 66.9 2022 7900 No es un Mazda2 VENCAMBIO RENAULT STEPWAY INTENS AUTOMÁTICO 2022 RECIBO KIA
30 119.9 2021 29000 No es un Mazda2 Vencambio mazda cx30 modelo 2021 grand touring lx refull 4×2
80 36.0 2011 117600 No es un Mazda2 Mazda 3 AllNew 2011. 2 litros. Full equipo.Los Alpes, Barran
105 22.0 2008 136 Salvamento HERMOSO MAZDA 2 EN VETA PARA USO PERSONAL (SALVAMENTO)Los Pi
106 23.5 2008 136000 Salvamento MAZDA 2 - 2008 SALVAMENTOMarsella Aloha Sector Sur, Kennedy2
113 22.0 2008 136 Duplicado Vendo hermoso Mazda 2 modelo 2008 (SALVAMENMarsella Aloha Se
189 29.9 1995 280000 No es un Mazda2 Mazda Estacas 4 x 2, buen estado gas gasolinaÁrea Urbana, Ch
203 28.0 2009 0 Km = 0 en carro usado Se vende Mazda 2Remanso, Villavicencio16 ago
255 32.0 2007 115000 No es un Mazda2 Vehículo Mazda 3 hatchback 115.000 Km 2 dueño, Excelente es
271 159.9 2022 0 No es un Mazda2 Mazda Cx5 Grand Touring 4*2 Plata | 2023Ciudad Campestre, Ca
273 169.8 2019 36026 No es un Mazda2 MAZDA CX-9 2019 cl 2.5 4x4 7 psj Blindaje 2Prados del Lim
280 6.8 2022 17200 Precio imposible Se vende hermoso Mazda 2 turing como nuevoLa Playa, Cali28 a
285 35.0 2007 148000 No es un Mazda2 Mazda 6 blindaje 2 plus modelo 2007 mecánico full equipo,co
312 107.5 2021 17909 No es un Mazda2 Venta Hermosa Camioneta Mazda Touring 2.o ATSan Pedro, CaliA 
mazda %>%
  filter(motivo == "Se conserva", km_en_miles) %>%
  transmute(id, modelo, km_reportado = kilometraje, km_corregido = kilometraje * 1000, anuncio) %>%
  kable(caption = "Kilometrajes escritos en miles (se multiplican por 1.000)") %>%
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)
Kilometrajes escritos en miles (se multiplican por 1.000)
id modelo km_reportado km_corregido anuncio
17 2011 183 183000 Mazda 2 1.5HM8La Mota, Medellín22 ago
168 2008 174 174000 SE VENDE Mazda 2La Ceja, Antioquia21 ago
180 2012 92 92000 hermoso carro Mazda 2 modelo 2012Aipe, Huilahace 6 días
185 2008 133 133000 MAZDA 2 MODELO 2008Yopal, Casanarehace 4 días
196 2019 54 54000 Mazda 2 prime sedanMira Rio II Etapa, Neiva02 sept
228 2013 116 116000 Vendo lindo Mazda 2 automáticoBavaria, Pereira29 ago
236 2010 114 114000 SE VENDE CARRO MAZDA 2 AUTOMÁTICOPaz Del Río, PereiraAyer
248 2019 51 51000 MAZDA 2 SEDAN GRAND TOURINGCentro, Bucaramanga04 sept
272 2011 99 99000 Mazda 2 mecanico 2011 Mecanico full equipo.Cañaverales - Los
298 2008 129 129000 Mazda 2 - 2008 automático en cueroEl Sena, Calihace 6 días
304 2019 35 35000 MAZDA 2 TOURING, modelo 2019. Mecánico, ÚNICO DUEÑOLili, Cal

Con esto construí la base de trabajo y las variables del análisis:

datos <- mazda %>%
  filter(motivo == "Se conserva") %>%
  mutate(
    kilometraje = ifelse(km_en_miles, kilometraje * 1000, kilometraje),
    modelo = ifelse(id == 299, 2009, modelo),
    precio_m = precio / 1e6,                 # precio en millones de pesos
    edad = 2022 - modelo,                    # años de antigüedad
    km_mil = kilometraje / 1000,             # miles de km
    transmision = factor(ifelse(grepl("^Autom", transmision), "Automática", "Mecánica"),
                         levels = c("Mecánica", "Automática")),
    generacion = factor(ifelse(modelo >= 2016, "Nueva (2016 en adelante)", "Anterior (2008 a 2015)"),
                        levels = c("Anterior (2008 a 2015)", "Nueva (2016 en adelante)")),
    region = case_when(
      departamento == "Antioquia" ~ "Antioquia",
      grepl("^Bogot", departamento) | departamento == "Cundinamarca" ~ "Bogotá y Cundinamarca",
      grepl("Valle", departamento) ~ "Valle del Cauca",
      departamento %in% c("Caldas", "Risaralda") | grepl("^Quind", departamento) ~ "Eje Cafetero",
      TRUE ~ "Otros departamentos"
    ),
    color_g = ifelse(color %in% c("Gris", "Blanco", "Rojo"), color, "Otro")
  )

n_final <- nrow(datos)
n_final
## [1] 299

Quedan 299 anuncios (se eliminaron 15). Además de las variables originales creé:

  • generacion: el Mazda2 cambió de generación en Colombia con el modelo 2016, que llegó al país en 2015 con la tecnología Skyactiv y un diseño nuevo (Autocosmos, 2015). Por eso separé los carros anteriores y posteriores a 2016.
  • region y color_g: agrupan departamentos y colores para no tener categorías con uno o dos anuncios.

Valores atípicos

# regla de Tukey sobre toda la base
q_precio <- quantile(datos$precio_m, c(0.25, 0.75))
lim_inf <- q_precio[1] - 1.5 * IQR(datos$precio_m)
lim_sup <- q_precio[2] + 1.5 * IQR(datos$precio_m)
n_atip_global <- sum(datos$precio_m < lim_inf | datos$precio_m > lim_sup)
n_atip_km <- sum(datos$km_mil > quantile(datos$km_mil, 0.75) + 1.5 * IQR(datos$km_mil))

# la misma regla, pero dentro de cada año del modelo
datos <- datos %>%
  group_by(modelo) %>%
  mutate(atipico_anio = precio_m < quantile(precio_m, 0.25) - 1.5 * IQR(precio_m) |
                        precio_m > quantile(precio_m, 0.75) + 1.5 * IQR(precio_m)) %>%
  ungroup()

ggplot(datos, aes(x = factor(modelo), y = precio_m, fill = generacion)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 2.5, outlier.fill = "white") +
  scale_fill_manual(values = c("steelblue", "darkorange")) +
  labs(title = "Precio anunciado según el año del modelo",
       x = "Año del modelo", y = "Precio (millones de pesos)", fill = "Generación") +
  theme(legend.position = "top")

datos %>%
  filter(atipico_anio) %>%
  select(id, modelo, precio_m, km_mil, transmision) %>%
  kable(digits = 1, caption = "Precios atípicos dentro de su año (regla de 1,5 veces el RIC)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Precios atípicos dentro de su año (regla de 1,5 veces el RIC)
id modelo precio_m km_mil transmision
10 2012 27.5 116.9 Mecánica
77 2011 22.0 102.0 Mecánica
78 2016 47.0 80.0 Mecánica
123 2017 28.9 12.3 Mecánica
179 2015 50.0 33.0 Mecánica
274 2017 72.0 81.0 Automática
295 2012 29.0 145.0 Mecánica
299 2009 25.8 129.0 Automática
306 2010 26.0 176.0 Mecánica 
n_atip_anio <- sum(datos$atipico_anio)
n_atip_mec <- sum(datos$atipico_anio & datos$transmision == "Mecánica")

Si se mira el precio en toda la base, no hay atípicos con la regla de 1,5 veces el rango intercuartílico (0 casos; el límite superior es 104.2 millones). Tampoco en el kilometraje (0 casos). Los carros más caros son sencillamente los más nuevos.

Como el precio depende tanto del año, tiene más sentido buscar atípicos dentro de cada año. Así aparecen 9 anuncios (los círculos blancos del gráfico), de los cuales 7 son mecánicos y varios tienen más kilómetros de lo normal para su año. Es decir, la mayoría tiene una razón para estar por debajo de los demás. El caso más raro es el id 123: un 2017 Touring con “12.345 km” (parece un número de relleno) anunciado a 28,9 millones, casi la mitad de lo que piden por otros 2017.

Decidí no eliminar ninguno de estos anuncios. No hay evidencia de que sean errores, y quitar datos solo porque se alejan del promedio haría que el modelo pareciera mejor de lo que es. Su efecto se revisa más adelante en la validación cruzada.

1.2 Análisis exploratorio

num <- datos %>% select(precio_m, edad, km_mil)

resumen <- data.frame(
  Media = sapply(num, mean),
  Mediana = sapply(num, median),
  DE = sapply(num, sd),
  Minimo = sapply(num, min),
  Q1 = sapply(num, quantile, 0.25),
  Q3 = sapply(num, quantile, 0.75),
  Maximo = sapply(num, max)
)
resumen$CV <- resumen$DE / resumen$Media
rownames(resumen) <- c("Precio (millones)", "Edad (años)", "Kilometraje (miles)")

kable(resumen, digits = 2, caption = "Estadísticas descriptivas de las variables cuantitativas") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Estadísticas descriptivas de las variables cuantitativas
Media Mediana DE Minimo Q1 Q3 Maximo CV
Precio (millones) 50.04 42 17.47 22 35.50 63 98.5 0.35
Edad (años) 7.06 7 4.06 0 4.00 10 14.0 0.57
Kilometraje (miles) 79.30 82 42.69 0 46.75 110 200.0 0.54 
bind_rows(
  datos %>% count(categoria = as.character(transmision)) %>% mutate(variable = "Transmisión"),
  datos %>% count(categoria = as.character(generacion)) %>% mutate(variable = "Generación"),
  datos %>% count(categoria = region) %>% mutate(variable = "Región"),
  datos %>% count(categoria = color_g) %>% mutate(variable = "Color")
) %>%
  group_by(variable) %>%
  mutate(porcentaje = round(100 * n / sum(n), 1)) %>%
  ungroup() %>%
  select(variable, categoria, n, porcentaje) %>%
  kable(caption = "Frecuencias de las variables categóricas") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Frecuencias de las variables categóricas
variable categoria n porcentaje
Transmisión Automática 164 54.8
Transmisión Mecánica 135 45.2
Generación Anterior (2008 a 2015) 167 55.9
Generación Nueva (2016 en adelante) 132 44.1
Región Antioquia 80 26.8
Región Bogotá y Cundinamarca 55 18.4
Región Eje Cafetero 43 14.4
Región Otros departamentos 71 23.7
Región Valle del Cauca 50 16.7
Color Blanco 76 25.4
Color Gris 101 33.8
Color Otro 79 26.4
Color Rojo 43 14.4 
datos %>%
  select(`Precio (millones)` = precio_m, `Edad (años)` = edad, `Kilometraje (miles)` = km_mil) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = valor)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white") +
  facet_wrap(~ variable, scales = "free") +
  labs(title = "Distribución de las variables cuantitativas", x = NULL, y = "Número de anuncios")

Algunas cosas que se ven en los datos:

  • El precio mediano es de 42 millones y la media de 50 millones. Que la media sea mayor indica una distribución sesgada a la derecha: muchos carros de 2011 a 2013 alrededor de 35 millones y una cola de carros recientes por encima de 80 millones. El coeficiente de variación es de 35 %, así que hay bastante variabilidad que explicar.
  • La edad mediana es de 7 años y el kilometraje mediano de 82 mil km. Hay carros desde 0 km (nuevos, publicados por concesionarios) hasta 200 mil km.
  • El 55 % de los anuncios son automáticos y el 44 % son de la generación nueva.
  • En el gráfico por año se ve que el precio baja con la antigüedad, pero entre los modelos 2015 y 2016 hay un salto que coincide con el cambio de generación.

1.3 Relación entre el precio y los predictores

Variables cuantitativas

datos %>%
  select(precio_m, generacion, `Edad (años)` = edad, `Kilometraje (miles de km)` = km_mil) %>%
  pivot_longer(c(`Edad (años)`, `Kilometraje (miles de km)`), names_to = "predictor", values_to = "valor") %>%
  ggplot(aes(x = valor, y = precio_m)) +
  geom_point(aes(color = generacion), alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "black", linetype = "dashed") +
  geom_smooth(method = "loess", formula = y ~ x, se = FALSE, color = "red") +
  facet_wrap(~ predictor, scales = "free_x") +
  scale_color_manual(values = c("steelblue", "darkorange")) +
  labs(title = "Precio frente a edad y kilometraje",
       subtitle = "Línea punteada: recta de regresión. Línea roja: suavizado loess",
       x = NULL, y = "Precio (millones de pesos)", color = "Generación") +
  theme(legend.position = "top")

kable(cor(num), digits = 3, caption = "Correlación de Pearson") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Correlación de Pearson
precio_m edad km_mil
precio_m 1.000 -0.945 -0.851
edad -0.945 1.000 0.878
km_mil -0.851 0.878 1.000 
kable(cor(num, method = "spearman"), digits = 3, caption = "Correlación de Spearman") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Correlación de Spearman
precio_m edad km_mil
precio_m 1.000 -0.959 -0.858
edad -0.959 1.000 0.880
km_mil -0.858 0.880 1.000 
cor.test(datos$precio_m, datos$edad)
## 
##  Pearson's product-moment correlation
## 
## data:  datos$precio_m and datos$edad
## t = -49.767, df = 297, p-value < 0.00000000000000022
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9559168 -0.9313433
## sample estimates:
##       cor 
## -0.944947
cor.test(datos$precio_m, datos$km_mil)
## 
##  Pearson's product-moment correlation
## 
## data:  datos$precio_m and datos$km_mil
## t = -27.878, df = 297, p-value < 0.00000000000000022
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8791970 -0.8158831
## sample estimates:
##        cor 
## -0.8505944
r_edad <- cor(datos$precio_m, datos$edad)
r_km <- cor(datos$precio_m, datos$km_mil)
r_edad_km <- cor(datos$edad, datos$km_mil)

La edad tiene una correlación de -0.945 con el precio, una relación negativa muy fuerte. Aun así, la curva loess muestra que la caída no es una línea recta: en pesos, el carro pierde más valor en los primeros años que en los últimos. Esto sugiere que más adelante convendría probar una transformación (por ejemplo, el logaritmo del precio).

El kilometraje también se relaciona negativamente con el precio (r = -0.851), pero está muy correlacionado con la edad (r = 0.878): los carros más viejos han rodado más. Por eso parte de la relación entre precio y kilometraje es en realidad el efecto de la edad. La correlación de Spearman da resultados parecidos, de modo que estas conclusiones no dependen de unos pocos valores extremos.

Variables categóricas

datos %>%
  transmute(precio_m,
            `Transmisión` = as.character(transmision),
            `Generación` = as.character(generacion),
            `Región` = region,
            Color = color_g) %>%
  pivot_longer(-precio_m, names_to = "variable", values_to = "categoria") %>%
  ggplot(aes(x = categoria, y = precio_m)) +
  geom_boxplot(fill = "lightblue") +
  facet_wrap(~ variable, scales = "free_x", ncol = 2) +
  labs(title = "Precio por grupos", x = NULL, y = "Precio (millones de pesos)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

datos %>%
  group_by(transmision) %>%
  summarise(n = n(), mediana = median(precio_m), media = mean(precio_m), de = sd(precio_m),
            edad_promedio = mean(edad)) %>%
  kable(digits = 1, caption = "Precio y edad según la transmisión") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Precio y edad según la transmisión
transmision n mediana media de edad_promedio
Mecánica 135 39 44.2 12.4 7.9
Automática 164 57 54.9 19.5 6.3 
prueba_t <- t.test(precio_m ~ transmision, data = datos)
prueba_w <- wilcox.test(precio_m ~ transmision, data = datos)
prueba_gen <- wilcox.test(precio_m ~ generacion, data = datos)
prueba_region <- kruskal.test(precio_m ~ region, data = datos)
prueba_color <- kruskal.test(precio_m ~ color_g, data = datos)

data.frame(
  Variable = c("Transmisión", "Transmisión", "Generación", "Región", "Color"),
  Prueba = c("t de Welch", "Wilcoxon", "Wilcoxon", "Kruskal-Wallis", "Kruskal-Wallis"),
  p_valor = c(prueba_t$p.value, prueba_w$p.value, prueba_gen$p.value,
              prueba_region$p.value, prueba_color$p.value)
) %>%
  mutate(p_valor = format.pval(p_valor, digits = 3, eps = 0.001)) %>%
  kable(caption = "Pruebas de diferencia entre grupos") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Pruebas de diferencia entre grupos
Variable Prueba p_valor
Transmisión t de Welch < 0.001
Transmisión Wilcoxon < 0.001
Generación Wilcoxon < 0.001
Región Kruskal-Wallis 0.67455
Color Kruskal-Wallis 0.00146 
datos %>%
  group_by(modelo) %>%
  summarise(prop_auto = mean(transmision == "Automática")) %>%
  ggplot(aes(x = modelo, y = prop_auto)) +
  geom_col(fill = "darkorange") +
  scale_x_continuous(breaks = 2008:2022) +
  scale_y_continuous(labels = function(x) paste0(100 * x, "%"), limits = c(0, 1)) +
  labs(title = "Proporción de carros automáticos por año del modelo", x = "Año del modelo", y = NULL)

Los automáticos tienen un precio mediano de 57 millones y los mecánicos de 39 millones, y la diferencia es significativa (p < 0.001 con Wilcoxon). Pero hay que tener cuidado al interpretarla: los automáticos son en promedio más nuevos (6.3 años frente a 7.9), y el último gráfico muestra que la proporción de automáticos crece en los modelos recientes. Una parte de la diferencia de precio se debe a la edad y no a la caja.

La generación nueva tiene un precio mediano 1.8 veces el de la anterior (p < 0.001), aunque también son carros más nuevos. Para la región (p = 0.675) y el color (p = 0.001) las diferencias son menores y tampoco tienen en cuenta la edad, así que no las tomo como efectos propios de esas variables.

1.4 Modelos de regresión lineal simple

Planteo tres modelos, cada uno con un predictor distinto:

  • S1: precio en función de la edad.
  • S2: precio en función del kilometraje (en miles de km).
  • S3: precio en función de la transmisión, que es una variable binaria. La categoría de referencia es Mecánica.
s1 <- lm(precio_m ~ edad, data = datos)
s2 <- lm(precio_m ~ km_mil, data = datos)
s3 <- lm(precio_m ~ transmision, data = datos)

summary(s1)
## 
## Call:
## lm(formula = precio_m ~ edad, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -29.5168  -3.8350  -0.1478   3.4487  19.7469 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 78.75309    0.66524  118.38 <0.0000000000000002 ***
## edad        -4.06725    0.08173  -49.77 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.725 on 297 degrees of freedom
## Multiple R-squared:  0.8929, Adjusted R-squared:  0.8926 
## F-statistic:  2477 on 1 and 297 DF,  p-value: < 0.00000000000000022
summary(s2)
## 
## Call:
## lm(formula = precio_m ~ km_mil, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.436  -6.485  -0.553   5.498  24.733 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 77.63175    1.12374   69.08 <0.0000000000000002 ***
## km_mil      -0.34798    0.01248  -27.88 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.199 on 297 degrees of freedom
## Multiple R-squared:  0.7235, Adjusted R-squared:  0.7226 
## F-statistic: 777.2 on 1 and 297 DF,  p-value: < 0.00000000000000022
summary(s3)
## 
## Call:
## lm(formula = precio_m ~ transmision, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.064 -12.769  -4.174  13.136  43.636 
## 
## Coefficients:
##                       Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)             44.174      1.434  30.806 < 0.0000000000000002 ***
## transmisionAutomática   10.690      1.936   5.521         0.0000000736 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.66 on 297 degrees of freedom
## Multiple R-squared:  0.09308,    Adjusted R-squared:  0.09002 
## F-statistic: 30.48 on 1 and 297 DF,  p-value: 0.0000000736
modelos <- list(s1, s2, s3)

tabla_simples <- data.frame(
  Modelo = c("S1: precio ~ edad", "S2: precio ~ km", "S3: precio ~ transmisión"),
  Intercepto = sapply(modelos, function(m) coef(m)[1]),
  Pendiente = sapply(modelos, function(m) coef(m)[2]),
  IC95_inf = sapply(modelos, function(m) confint(m)[2, 1]),
  IC95_sup = sapply(modelos, function(m) confint(m)[2, 2]),
  p_valor = sapply(modelos, function(m) format.pval(summary(m)$coefficients[2, 4], digits = 3, eps = 0.001)),
  R2 = sapply(modelos, function(m) summary(m)$r.squared),
  Error_residual = sapply(modelos, function(m) summary(m)$sigma)
)

kable(tabla_simples, digits = 3, caption = "Resumen de los modelos simples (precio en millones)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Resumen de los modelos simples (precio en millones)
Modelo Intercepto Pendiente IC95_inf IC95_sup p_valor R2 Error_residual
edad S1: precio ~ edad 78.753 -4.067 -4.228 -3.906 <0.001 0.893 5.725
km_mil S2: precio ~ km 77.632 -0.348 -0.373 -0.323 <0.001 0.724 9.199
transmisionAutomática S3: precio ~ transmisión 44.174 10.690 6.879 14.500 <0.001 0.093 16.661

1.5 Interpretación de los modelos

S1 (edad). El intercepto, 78.8 millones, es el precio estimado de un carro con edad 0, es decir, un modelo 2022. La pendiente indica que cada año adicional de antigüedad se asocia, en promedio, con un precio anunciado 4.07 millones menor. La edad sola explica el 89.3 % de la variación del precio y el coeficiente es claramente significativo (p < 0.001).

S2 (kilometraje). El intercepto (77.6 millones) sería el precio de un carro sin recorrido. Por cada 10.000 km adicionales el precio anunciado es, en promedio, 3.48 millones menor (p < 0.001). Su R² es de 72.4 %, bastante menor que el de la edad.

S3 (transmisión). El intercepto (44.2 millones) es el precio promedio de los carros mecánicos, que son la categoría de referencia. El coeficiente indica que los automáticos se anuncian en promedio 10.69 millones más caros (p < 0.001). La diferencia es significativa, pero el R² es apenas 9.3 %: la transmisión por sí sola explica muy poco del precio.

Los tres coeficientes son significativos al 5 %, pero eso no quiere decir que los tres modelos sean igual de útiles. La significancia solo indica que la relación probablemente no es cero. Para predecir importa más cuánto se equivoca el modelo, que es lo que se evalúa en la sección 1.6.

Asociación no es causalidad. Estos coeficientes describen cómo se relaciona el precio con cada variable en estos anuncios, no cuánto cambiaría el precio de un carro si cambiara solo esa característica. El kilometraje es un buen ejemplo: en S2 cada 10.000 km se asocian con 3.48 millones menos, pero si en el mismo modelo se incluye también la edad, ese valor baja a 0.38 millones. La mayor parte de lo que S2 le atribuye al kilometraje es en realidad la edad del carro. Lo mismo pasa con la transmisión, porque los automáticos tienden a ser más nuevos. Además hay variables que no están en la base (estado del carro, versión, número de dueños, accidentes) y que también afectan el precio. Con datos observacionales como estos solo puedo hablar de asociaciones.

Antes de comparar los modelos reviso los residuos de S1, que es el de mayor R²:

data.frame(edad = datos$edad, residuo = residuals(s1)) %>%
  ggplot(aes(x = factor(edad), y = residuo)) +
  geom_hline(yintercept = 0, color = "gray40") +
  geom_boxplot(fill = "lightblue") +
  labs(title = "Residuos de S1 según la edad del carro",
       x = "Edad (años)", y = "Residuo (millones)")

Los residuos no están centrados en cero para todas las edades. El modelo subestima los carros del año y los más viejos (12 a 14 años) y sobreestima los de 7 a 10 años (modelos 2012 a 2015), lo que confirma que la relación es curva. Entre los 6 y los 7 años hay además un salto brusco, justo donde cambia la generación del Mazda2. La dispersión también es mayor en los carros de 0 a 5 años que en los viejos. Son problemas que se trabajarán en el Punto 2 con modelos múltiples y una transformación de la respuesta.

1.6 Validación cruzada y elección del mejor modelo simple

Para comparar los modelos por su capacidad de predecir carros que no se usaron para estimarlos, uso validación cruzada de 10 particiones repetida 10 veces. En cada repetición los datos se dividen al azar en 10 grupos; el modelo se estima con 9 y se evalúa con el grupo restante. Esto da 100 estimaciones del error por modelo. Algunas condiciones:

  • Los tres modelos se evalúan con las mismas particiones.
  • La métrica es el error absoluto medio (MAE), en millones de pesos.
  • Como referencia uso un modelo que siempre predice la mediana del precio en los datos de entrenamiento. Si un modelo no mejora esa referencia, no sirve.
formulas <- list(
  "S1: edad" = precio_m ~ edad,
  "S2: kilometraje" = precio_m ~ km_mil,
  "S3: transmisión" = precio_m ~ transmision
)

correr_cv <- function(semilla, R = 10, K = 10) {
  set.seed(semilla)
  n <- nrow(datos)
  resultados <- data.frame()
  for (r in 1:R) {
    grupos <- sample(rep(1:K, length.out = n))
    for (k in 1:K) {
      train <- datos[grupos != k, ]
      test <- datos[grupos == k, ]
      fila <- data.frame(rep = r, fold = k,
                         Mediana = mean(abs(test$precio_m - median(train$precio_m))))
      for (nombre in names(formulas)) {
        m <- lm(formulas[[nombre]], data = train)
        fila[[nombre]] <- mean(abs(test$precio_m - predict(m, newdata = test)))
      }
      resultados <- rbind(resultados, fila)
    }
  }
  resultados
}

cv <- correr_cv(semilla = 2026)

mae_cv <- cv %>%
  select(-rep, -fold) %>%
  pivot_longer(everything(), names_to = "Modelo", values_to = "MAE") %>%
  group_by(Modelo) %>%
  summarise(MAE_promedio = mean(MAE), DE = sd(MAE), Minimo = min(MAE), Maximo = max(MAE)) %>%
  mutate(Error_relativo = MAE_promedio / median(datos$precio_m),
         Mejora_vs_mediana = 1 - MAE_promedio / MAE_promedio[Modelo == "Mediana"]) %>%
  arrange(MAE_promedio)

kable(mae_cv, digits = 3,
      caption = "MAE en validación cruzada 10 x 10 (millones de pesos). Error relativo = MAE / precio mediano") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
MAE en validación cruzada 10 x 10 (millones de pesos). Error relativo = MAE / precio mediano
Modelo MAE_promedio DE Minimo Maximo Error_relativo Mejora_vs_mediana
S1: edad 4.390 0.590 3.182 6.098 0.105 0.702
S2: kilometraje 7.208 1.235 4.813 10.575 0.172 0.510
S3: transmisión 14.400 1.475 10.245 18.416 0.343 0.021
Mediana 14.714 1.761 11.057 20.573 0.350 0.000 
cv %>%
  select(-rep, -fold) %>%
  pivot_longer(everything(), names_to = "Modelo", values_to = "MAE") %>%
  mutate(Modelo = factor(Modelo, levels = rev(mae_cv$Modelo))) %>%
  ggplot(aes(x = MAE, y = Modelo)) +
  geom_boxplot(fill = "lightblue") +
  labs(title = "MAE en cada una de las 100 particiones",
       x = "MAE (millones de pesos)", y = NULL)

Predecir siempre la mediana tiene un error promedio de 14.71 millones. S1 lo reduce a 4.39 millones, una mejora de 70 %, y su error equivale al 10.5 % del precio mediano. S2 queda en 7.21 millones. S3 apenas mejora la referencia (14.4 millones) y solo le gana a la mediana en el 63 % de las particiones.

Comparando partición por partición, S1 tiene menor error que S2 en el 99 % de los casos, así que la diferencia entre ambos no depende de cómo se dividieron los datos.

Verificaciones adicionales

Para comprobar que el resultado no depende de la semilla, repetí la validación cruzada con cinco semillas distintas (5 repeticiones de 10 particiones cada una):

semillas <- c(1, 42, 123, 777, 2024)
sens <- t(sapply(semillas, function(s) colMeans(correr_cv(s, R = 5)[, -(1:2)])))
sens <- data.frame(Semilla = semillas, sens, check.names = FALSE)

kable(sens, digits = 3, caption = "MAE promedio con distintas semillas (millones de pesos)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
MAE promedio con distintas semillas (millones de pesos)
Semilla Mediana S1: edad S2: kilometraje S3: transmisión
1 14.712 4.385 7.218 14.436
42 14.711 4.392 7.215 14.423
123 14.713 4.393 7.211 14.402
777 14.698 4.394 7.203 14.408
2024 14.706 4.396 7.217 14.398 
orden_igual <- all(sens$`S1: edad` < sens$`S2: kilometraje` &
                   sens$`S2: kilometraje` < sens$`S3: transmisión` &
                   sens$`S3: transmisión` < sens$Mediana)

Con todas las semillas el orden es el mismo: S1, luego S2, luego S3 y por último la mediana. Los valores cambian solo en el segundo decimal.

También calculé el error con validación cruzada dejando uno por fuera (LOOCV). En regresión lineal no hace falta reajustar el modelo n veces: el residuo de cada observación cuando se deja por fuera es ei / (1 - hii), donde hii es su apalancamiento.

mae_loo <- function(m) mean(abs(residuals(m) / (1 - hatvalues(m))))
y <- datos$precio_m
mae_loo_mediana <- mean(sapply(1:length(y), function(i) abs(y[i] - median(y[-i]))))

data.frame(
  Modelo = c("S1: edad", "S2: kilometraje", "S3: transmisión", "Mediana"),
  MAE_LOOCV = c(mae_loo(s1), mae_loo(s2), mae_loo(s3), mae_loo_mediana),
  MAE_CV_10x10 = c(mae_s1, mae_s2, mae_s3, mae_med)
) %>%
  kable(digits = 3, caption = "Comparación entre LOOCV y validación cruzada 10 x 10") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Comparación entre LOOCV y validación cruzada 10 x 10
Modelo MAE_LOOCV MAE_CV_10x10
S1: edad 4.391 4.390
S2: kilometraje 7.212 7.208
S3: transmisión 14.403 14.400
Mediana 14.705 14.714

Los dos métodos dan prácticamente el mismo resultado.

Elección del mejor modelo simple

El mejor modelo simple es S1 (precio en función de la edad):

  • Precisión: es el que menos se equivoca, con un MAE de 4.39 millones, y eso se mantiene con distintas semillas y con LOOCV.
  • Coherencia: el año del modelo es lo primero que mira cualquier comprador o vendedor de carros usados, y el signo negativo tiene sentido.
  • Interpretabilidad: el coeficiente se lee directamente como millones de pesos por año de antigüedad.

Su principal limitación es que supone que el carro pierde la misma cantidad de pesos cada año, cuando los datos muestran una caída más fuerte en los primeros años y un salto entre las dos generaciones. Además, un error promedio de 4.4 millones todavía es alto para fijar el precio de un carro concreto. Por eso en el Punto 2 se combinan varias variables y se prueba una transformación del precio.

Conclusiones del punto 1

  • La base necesitó una limpieza cuidadosa: 15 de los 314 anuncios no eran comparables (otros modelos de Mazda, un carro de salvamento repetido, precios o kilometrajes imposibles), y 12 datos se corrigieron con la información del mismo anuncio.
  • La edad es la variable más relacionada con el precio (r = -0.94). El kilometraje también lo está, pero en gran parte porque los carros viejos tienen más kilómetros. La transmisión muestra una diferencia de precio que se explica en parte porque los automáticos son más nuevos.
  • De los tres modelos simples, el de la edad es el mejor en validación cruzada (MAE de 4.39 millones frente a 14.71 de la referencia), y el resultado se mantiene con otras semillas y con LOOCV.
  • Los residuos de S1 muestran una relación curva y un cambio entre generaciones, lo que indica que el modelo simple se queda corto y justifica los modelos del Punto 2.
  • Todo lo anterior se refiere al precio anunciado en 2022, no al precio real de venta.