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

library(readxl)
library(MASS)       # stepAIC (se carga antes de dplyr para no tapar select)
library(car)        # vif
library(lmtest)     # bptest, dwtest
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.

Punto 1. Exploración y regresión lineal simple

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:

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.
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.
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).
Un carro 2009 tiene 0 km, es decir, no se reportó el kilometraje.
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).
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 las verificaciones adicionales del Punto 2 (sección 2.6).

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:

Todos los modelos del trabajo (los simples de este punto y los múltiples del Punto 2) se evalúan con las mismas particiones, que se crean una sola vez.
La métrica principal es el error absoluto medio (MAE), en millones de pesos. También reporto el RMSE, el MAPE (error en porcentaje del precio) y la cobertura del intervalo de predicción del 95 %, es decir, qué porcentaje de los precios reales queda dentro del intervalo.
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.

crear_particiones <- function(semilla, R = 10, K = 10) {
  set.seed(semilla)
  sapply(1:R, function(r) sample(rep(1:K, length.out = nrow(datos))))
}

# estima el modelo con 9 grupos y predice el grupo restante, para todas las particiones
cv_modelo <- function(formula, particiones, log = FALSE, paso_aic = FALSE) {
  f <- as.formula(formula)
  salida <- list()
  for (r in 1:ncol(particiones)) {
    for (k in 1:max(particiones)) {
      train <- datos[particiones[, r] != k, ]
      test <- datos[particiones[, r] == k, ]
      m <- do.call("lm", list(formula = f, data = train))
      if (paso_aic) m <- stepAIC(m, direction = "both", trace = FALSE)  # selección dentro del entrenamiento
      p <- predict(m, newdata = test, interval = "prediction", level = 0.95)
      if (log) p <- exp(p)   # regreso a millones de pesos
      salida[[length(salida) + 1]] <- data.frame(
        rep = r, fold = k, id = test$id, real = test$precio_m,
        pred = p[, "fit"], inf = p[, "lwr"], sup = p[, "upr"],
        variables = paste(attr(terms(m), "term.labels"), collapse = " + "))
    }
  }
  do.call(rbind, salida)
}

cv_mediana <- function(particiones) {
  salida <- list()
  for (r in 1:ncol(particiones)) {
    for (k in 1:max(particiones)) {
      train <- datos[particiones[, r] != k, ]
      test <- datos[particiones[, r] == k, ]
      salida[[length(salida) + 1]] <- data.frame(
        rep = r, fold = k, id = test$id, real = test$precio_m,
        pred = median(train$precio_m), inf = NA, sup = NA, variables = "mediana")
    }
  }
  do.call(rbind, salida)
}

mae_por_particion <- function(res) {
  res %>%
    group_by(rep, fold) %>%
    summarise(mae = mean(abs(real - pred)), .groups = "drop") %>%
    pull(mae)
}

resumen_cv <- function(lista) {
  data.frame(
    Modelo = names(lista),
    MAE = sapply(lista, function(x) mean(mae_por_particion(x))),
    DE_MAE = sapply(lista, function(x) sd(mae_por_particion(x))),
    RMSE = sapply(lista, function(x) sqrt(mean((x$real - x$pred)^2))),
    MAPE = sapply(lista, function(x) 100 * mean(abs(x$real - x$pred) / x$real)),
    Cobertura_IP95 = sapply(lista, function(x) 100 * mean(x$real >= x$inf & x$real <= x$sup)),
    row.names = NULL
  ) %>%
    mutate(Mejora_vs_mediana = 100 * (1 - MAE / MAE[Modelo == "Mediana"])) %>%
    arrange(MAE)
}

particiones <- crear_particiones(2026)

cv_res <- list(
  "Mediana" = cv_mediana(particiones),
  "S1: edad" = cv_modelo("precio_m ~ edad", particiones),
  "S2: kilometraje" = cv_modelo("precio_m ~ km_mil", particiones),
  "S3: transmisión" = cv_modelo("precio_m ~ transmision", particiones)
)

tabla_cv_simples <- resumen_cv(cv_res)
kable(tabla_cv_simples, digits = 2,
      caption = "Validación cruzada 10 x 10 de los modelos simples (errores en millones de pesos; MAPE, cobertura y mejora en %)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Validación cruzada 10 x 10 de los modelos simples (errores en millones de pesos; MAPE, cobertura y mejora en %)
Modelo	MAE	DE_MAE	RMSE	MAPE	Cobertura_IP95	Mejora_vs_mediana
S1: edad	4.39	0.59	5.75	9.66	96.32	70.17
S2: kilometraje	7.21	1.23	9.24	16.18	94.98	51.01
S3: transmisión	14.40	1.48	16.71	31.32	96.89	2.13
Mediana	14.71	1.76	19.12	26.79	NA	0.00

mae_largo <- bind_rows(lapply(names(cv_res), function(nm) {
  data.frame(Modelo = nm, MAE = mae_por_particion(cv_res[[nm]]))
}))
mae_largo$Modelo <- factor(mae_largo$Modelo, levels = rev(tabla_cv_simples$Modelo))

ggplot(mae_largo, 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.

Verificación con otras semillas

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_simples <- data.frame()
for (s in semillas) {
  part_s <- crear_particiones(s, R = 5)
  sens_simples <- rbind(sens_simples, data.frame(
    Semilla = s,
    Mediana = mean(mae_por_particion(cv_mediana(part_s))),
    S1 = mean(mae_por_particion(cv_modelo("precio_m ~ edad", part_s))),
    S2 = mean(mae_por_particion(cv_modelo("precio_m ~ km_mil", part_s))),
    S3 = mean(mae_por_particion(cv_modelo("precio_m ~ transmision", part_s)))
  ))
}

kable(sens_simples, 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	S2	S3
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_simples$S1 < sens_simples$S2 & sens_simples$S2 < sens_simples$S3 &
                   sens_simples$S3 < sens_simples$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.

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.
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.
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.

Punto 2. Regresión lineal múltiple

2.1 Modelos propuestos

Propongo dos modelos múltiples:

M1 (básico): precio en función de la edad, el kilometraje y la transmisión. Son las tres características que cualquier comprador pregunta primero, y las tres mostraron relación con el precio en el Punto 1. Incluye una categórica binaria (transmisión, con referencia Mecánica).
M2 (ampliado): M1 más la generación, la región y el color.
- La generación recoge el cambio de producto de 2016, que la edad por sí sola no captura (es el salto que se vio en los residuos de S1). Referencia: generación anterior (2008 a 2015).
- La región se incluye porque el mercado de carros usados puede tener precios distintos según la ciudad. Referencia: Bogotá y Cundinamarca.
- El color entra como candidato, porque algunos colores pueden ser más buscados que otros. Referencia: Gris, el más frecuente.

datos <- datos %>%
  mutate(region = relevel(factor(region), ref = "Bogotá y Cundinamarca"),
         color_g = relevel(factor(color_g), ref = "Gris"))

levels(datos$transmision)

## [1] "Mecánica"   "Automática"

levels(datos$generacion)

## [1] "Anterior (2008 a 2015)"   "Nueva (2016 en adelante)"

levels(datos$region)

## [1] "Bogotá y Cundinamarca" "Antioquia"             "Eje Cafetero"         
## [4] "Otros departamentos"   "Valle del Cauca"

levels(datos$color_g)

## [1] "Gris"   "Blanco" "Otro"   "Rojo"

m1 <- lm(precio_m ~ edad + km_mil + transmision, data = datos)
m2 <- lm(precio_m ~ edad + km_mil + transmision + generacion + region + color_g, data = datos)

tabla_coef <- function(m) {
  s <- summary(m)$coefficients
  ic <- confint(m)
  data.frame(Termino = rownames(s), Estimacion = s[, 1], Error_est = s[, 2], t = s[, 3],
             p_valor = format.pval(s[, 4], digits = 3, eps = 0.001),
             IC95_inf = ic[, 1], IC95_sup = ic[, 2], row.names = NULL)
}

kable(tabla_coef(m1), digits = 3, caption = "M1: coeficientes (precio en millones)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

M1: coeficientes (precio en millones)
Termino	Estimacion	Error_est	t	p_valor	IC95_inf	IC95_sup
(Intercept)	76.106	0.786	96.821	<0.001	74.560	77.653
edad	-3.617	0.157	-23.009	<0.001	-3.927	-3.308
km_mil	-0.037	0.015	-2.507	0.0127	-0.067	-0.008
transmisionAutomática	4.428	0.623	7.108	<0.001	3.202	5.654

kable(tabla_coef(m2), digits = 3, caption = "M2: coeficientes (precio en millones)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

M2: coeficientes (precio en millones)
Termino	Estimacion	Error_est	t	p_valor	IC95_inf	IC95_sup
(Intercept)	66.107	1.697	38.956	< 0.001	62.767	69.447
edad	-2.758	0.191	-14.441	< 0.001	-3.134	-2.382
km_mil	-0.038	0.014	-2.719	0.00694	-0.065	-0.010
transmisionAutomática	3.868	0.591	6.548	< 0.001	2.705	5.030
generacionNueva (2016 en adelante)	8.356	1.179	7.090	< 0.001	6.037	10.676
regionAntioquia	2.540	0.875	2.903	0.00398	0.818	4.262
regionEje Cafetero	2.167	1.011	2.144	0.03287	0.178	4.155
regionOtros departamentos	0.524	0.890	0.588	0.55696	-1.229	2.276
regionValle del Cauca	0.681	0.963	0.708	0.47979	-1.214	2.576
color_gBlanco	-0.654	0.756	-0.865	0.38753	-2.142	0.834
color_gOtro	-1.068	0.751	-1.422	0.15613	-2.547	0.410
color_gRojo	-1.309	0.915	-1.431	0.15365	-3.109	0.492

2.2 Interpretación y selección paso a paso con AIC

M1. Manteniendo constantes el kilometraje y la transmisión, cada año adicional de antigüedad se asocia con un precio 3.62 millones menor. Lo más interesante es lo que pasa con el kilometraje: en S2 cada 10.000 km se asociaban con 3.48 millones menos, pero en M1, a igual edad, solo con 0.37 millones menos (p = 0.013). Algo parecido ocurre con la transmisión: la diferencia a favor de los automáticos baja de 10.69 millones en S3 a 4.43 millones en M1. Esto confirma lo que se veía en el Punto 1: buena parte de lo que los modelos simples le atribuían al kilometraje y a la transmisión era en realidad el efecto de la edad. El R² ajustado de M1 es 0.909.

M2. A igual edad, kilometraje y transmisión, un carro de la generación nueva se anuncia 8.36 millones más caro (p < 0.001). Frente a Bogotá y Cundinamarca, los carros de Antioquia se anuncian 2.54 millones más caros (p = 0.004). Ningún color se diferencia significativamente del gris. El R² ajustado sube a 0.924. El R² ajustado se usa en lugar del R² normal porque penaliza agregar variables que no aportan.

Ahora aplico la selección paso a paso con AIC sobre M2, en ambas direcciones (quitar y agregar variables):

m2_aic <- stepAIC(m2, direction = "both", trace = FALSE)
m2_aic$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## precio_m ~ edad + km_mil + transmision + generacion + region + 
##     color_g
## 
## Final Model:
## precio_m ~ edad + km_mil + transmision + generacion + region
## 
## 
##        Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                             287   6683.419 952.9755
## 2 - color_g  3 69.78065       290   6753.200 950.0811

kable(tabla_coef(m2_aic), digits = 3, caption = "M2-AIC: modelo elegido por el AIC") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

M2-AIC: modelo elegido por el AIC
Termino	Estimacion	Error_est	t	p_valor	IC95_inf	IC95_sup
(Intercept)	65.649	1.665	39.431	< 0.001	62.372	68.926
edad	-2.758	0.190	-14.481	< 0.001	-3.133	-2.383
km_mil	-0.040	0.014	-2.892	0.00412	-0.067	-0.013
transmisionAutomática	3.951	0.587	6.736	< 0.001	2.796	5.105
generacionNueva (2016 en adelante)	8.165	1.164	7.014	< 0.001	5.874	10.456
regionAntioquia	2.726	0.863	3.161	0.00174	1.029	4.424
regionEje Cafetero	2.071	1.005	2.061	0.04020	0.093	4.048
regionOtros departamentos	0.484	0.877	0.552	0.58132	-1.242	2.210
regionValle del Cauca	0.695	0.954	0.729	0.46686	-1.183	2.573

comparacion_lineal <- data.frame(
  Modelo = c("M1", "M2", "M2-AIC"),
  Parametros = sapply(list(m1, m2, m2_aic), function(m) length(coef(m))),
  R2 = sapply(list(m1, m2, m2_aic), function(m) summary(m)$r.squared),
  R2_ajustado = sapply(list(m1, m2, m2_aic), function(m) summary(m)$adj.r.squared),
  AIC = sapply(list(m1, m2, m2_aic), AIC),
  Error_residual = sapply(list(m1, m2, m2_aic), sigma)
)
kable(comparacion_lineal, digits = 3, caption = "Comparación de los modelos lineales (mismas observaciones y misma escala)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Comparación de los modelos lineales (mismas observaciones y misma escala)
Modelo	Parametros	R2	R2_ajustado	AIC	Error_residual
M1	4	0.910	0.909	1847.083	5.258
M2	12	0.926	0.924	1803.501	4.826
M2-AIC	9	0.926	0.924	1800.606	4.826

nombres_var <- c(edad = "edad", km_mil = "kilometraje", transmision = "transmisión",
                 generacion = "generación", region = "región", color_g = "color")
eliminadas <- setdiff(attr(terms(m2), "term.labels"), attr(terms(m2_aic), "term.labels"))
texto_eliminadas <- ifelse(length(eliminadas) == 0, "ninguna variable",
                           paste(nombres_var[eliminadas], collapse = " y "))

El procedimiento elimina color y el AIC baja de 1803.5 a 1800.6. Los coeficientes de las demás variables prácticamente no cambian, lo que indica que la variable eliminada no estaba alterando las otras relaciones. El AIC mantiene la región aunque no todas sus categorías sean significativas, porque su criterio es distinto al de los valores p: conserva una variable si mejora el ajuste lo suficiente para compensar los parámetros que agrega. En la validación cruzada esta selección se repite dentro de cada partición de entrenamiento, para que la evaluación no use información de los datos de prueba.

2.3 Supuestos del modelo (M2-AIC)

par(mfrow = c(2, 2))
plot(m2_aic)

par(mfrow = c(1, 1))

bp_lin <- bptest(m2_aic)
sw_lin <- shapiro.test(residuals(m2_aic))
dw_lin <- dwtest(m2_aic)   # en el orden en que se extrajeron los anuncios

tabla_pruebas <- function(bp, sw, dw) {
  data.frame(
    Supuesto = c("Varianza constante", "Normalidad de los errores", "Independencia"),
    Prueba = c("Breusch-Pagan", "Shapiro-Wilk", "Durbin-Watson"),
    Estadistico = c(bp$statistic, sw$statistic, dw$statistic),
    p_valor = format.pval(c(bp$p.value, sw$p.value, dw$p.value), digits = 3, eps = 0.001),
    Se_rechaza_al_5 = ifelse(c(bp$p.value, sw$p.value, dw$p.value) < 0.05, "Sí", "No")
  )
}
kable(tabla_pruebas(bp_lin, sw_lin, dw_lin), digits = 3, caption = "Pruebas de supuestos para M2-AIC") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Pruebas de supuestos para M2-AIC
	Supuesto	Prueba	Estadistico	p_valor	Se_rechaza_al_5
BP	Varianza constante	Breusch-Pagan	27.431	<0.001	Sí
W	Normalidad de los errores	Shapiro-Wilk	0.929	<0.001	Sí
DW	Independencia	Durbin-Watson	1.960	0.286	No

vif(m2_aic)

##                 GVIF Df GVIF^(1/(2*Df))
## edad        7.642023  1        2.764421
## km_mil      4.461156  1        2.112145
## transmision 1.093836  1        1.045866
## generacion  4.289332  1        2.071070
## region      1.154382  4        1.018108

vif_lin <- vif(m2_aic)

cook_lin <- cooks.distance(m2_aic)
data.frame(id = datos$id, cook = cook_lin) %>%
  ggplot(aes(x = id, y = cook)) +
  geom_col(fill = "steelblue") +
  geom_hline(yintercept = 4 / n_final, linetype = "dashed", color = "red") +
  labs(title = "Distancia de Cook de M2-AIC", subtitle = "Línea roja: 4/n",
       x = "id del anuncio", y = "Distancia de Cook")

Linealidad. En el gráfico de residuos contra valores ajustados (el primero de los cuatro) los residuos forman una curva en U: el modelo subestima los carros más baratos y los más caros y sobreestima los del medio. Si la relación fuera lineal, los puntos estarían repartidos parejos alrededor de cero. La consecuencia es que el modelo se equivoca de forma sistemática según el rango de precio.
Varianza constante. La dispersión de los residuos crece con el valor ajustado, y la prueba de Breusch-Pagan rechaza la homocedasticidad (p < 0.001). Es lógico: equivocarse un 5 % en un carro de 90 millones es más plata que en uno de 30. La consecuencia es que los errores estándar y los intervalos de predicción del modelo no son confiables, porque suponen el mismo ancho para todos los carros.
Normalidad. Shapiro-Wilk rechaza la normalidad (p < 0.001) y el gráfico Q-Q muestra colas más pesadas que las de una normal. Con 299 datos los coeficientes siguen siendo aproximadamente normales, pero los intervalos de predicción sí dependen de este supuesto.
Multicolinealidad. El VIF más alto es el de la edad (7.64), porque la edad está relacionada con el kilometraje y con la generación. Es una colinealidad moderada (menor a 10). Hace que los coeficientes de edad y kilometraje sean menos precisos, pero no daña las predicciones.
Independencia. El estadístico de Durbin-Watson es 1.96, cercano a 2, así que no hay autocorrelación en el orden en que se extrajeron los anuncios. No hay una dimensión de tiempo, porque todos son de la misma época. El riesgo que sí existe es que un mismo concesionario publique varios carros parecidos.
Observaciones influyentes. 13 anuncios superan el umbral de Cook de 4/n, pero ninguno tiene una distancia mayor a 1 (la máxima es 0.21). Se revisan en detalle en el modelo final.

En resumen, el modelo lineal tiene problemas de forma funcional y de varianza que se pueden corregir con una transformación.

2.4 Transformación de la respuesta

Hay dos razones para transformar el precio con logaritmo:

Un carro no pierde la misma cantidad de pesos cada año, sino más o menos el mismo porcentaje de su valor. Ese tipo de relación se representa mejor en escala logarítmica.
Los diagnósticos de M2-AIC muestran curvatura y varianza creciente, dos problemas que el logaritmo suele corregir.

Además, más adelante compruebo con la validación cruzada que la transformación realmente reduce el error en pesos.

Estimé los mismos modelos con log(precio) y apliqué otra vez la selección por AIC:

m1_log <- lm(log(precio_m) ~ edad + km_mil + transmision, data = datos)
m2_log <- lm(log(precio_m) ~ edad + km_mil + transmision + generacion + region + color_g, data = datos)
m2_log_aic <- stepAIC(m2_log, direction = "both", trace = FALSE)
m3_log <- lm(log(precio_m) ~ edad + km_mil + transmision + generacion, data = datos)

attr(terms(m2_log_aic), "term.labels")

## [1] "edad"        "km_mil"      "transmision" "generacion"  "region"

tabla_coef_log <- function(m) {
  t <- tabla_coef(m)
  t$Efecto_pct <- ifelse(t$Termino == "(Intercept)", NA, 100 * (exp(t$Estimacion) - 1))
  t
}
kable(tabla_coef_log(m3_log), digits = 4,
      caption = "M3-log: coeficientes en escala logarítmica y efecto en % sobre el precio") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

M3-log: coeficientes en escala logarítmica y efecto en % sobre el precio
Termino	Estimacion	Error_est	t	p_valor	IC95_inf	IC95_sup	Efecto_pct
(Intercept)	4.1850	0.0260	160.7707	<0.001	4.1337	4.2362	NA
edad	-0.0562	0.0032	-17.4084	<0.001	-0.0626	-0.0499	-5.4675
km_mil	-0.0005	0.0002	-2.0779	0.0386	-0.0010	0.0000	-0.0489
transmisionAutomática	0.0593	0.0099	6.0059	<0.001	0.0399	0.0787	6.1085
generacionNueva (2016 en adelante)	0.1655	0.0195	8.4740	<0.001	0.1271	0.2039	17.9991

En escala logarítmica, el AIC conserva estas variables: edad, kilometraje, transmisión, generación, región. Es decir, elimina el color, igual que antes, pero mantiene la región. También estimé M3-log, que es igual pero sin la región. Es más sencillo de usar, porque solo pide datos que cualquier dueño conoce de su carro, y no hay una razón clara para que un mismo Mazda2 valga distinto según la ciudad. En la sección 2.5 se compara si la región realmente mejora la predicción.

Interpretación de M3-log. Como la respuesta está en logaritmo, cada coeficiente se convierte en un cambio porcentual con la fórmula 100 x (exp(β) - 1):

Edad: cada año adicional se asocia con un precio 5.5 % menor, manteniendo lo demás constante. Es la depreciación anual promedio.
Kilometraje: cada 10.000 km adicionales se asocian con un precio 0.49 % menor (p = 0.039). Una vez se conoce la edad, el kilometraje mueve poco el precio anunciado.
Transmisión automática: se asocia con un precio 6.1 % mayor que el de un carro mecánico parecido.
Generación nueva: se asocia con un precio 18 % mayor que el de un carro de la generación anterior con la misma edad, kilometraje y transmisión.
El R² ajustado es 0.939, pero no se puede comparar con el de los modelos lineales, porque mide qué tanto se explica la variación del logaritmo del precio y no la del precio.

Sobre la comparación con el modelo lineal

Como indica el enunciado, no comparo el R² ni el AIC del modelo logarítmico con los del modelo lineal, porque la respuesta está en otra escala (logaritmo del precio frente a precio). La comparación justa entre los dos es el MAE en pesos, que se calcula más abajo con validación cruzada.

Supuestos del modelo transformado

corte <- which(sort(unique(datos$modelo)) == 2016) - 0.5
bind_rows(
  data.frame(modelo = datos$modelo, residuo = residuals(m1_log), Modelo = "M1-log (sin generación)"),
  data.frame(modelo = datos$modelo, residuo = residuals(m3_log), Modelo = "M3-log (con generación)")
) %>%
  ggplot(aes(x = factor(modelo), y = residuo)) +
  geom_hline(yintercept = 0, color = "gray40") +
  geom_vline(xintercept = corte, linetype = "dashed", color = "darkorange") +
  geom_boxplot(fill = "lightblue") +
  facet_wrap(~ Modelo, ncol = 1) +
  labs(title = "Residuos por año del modelo",
       subtitle = "Línea naranja: llegada de la generación nueva (modelo 2016)",
       x = "Año del modelo", y = "Residuo (escala log)")

Este gráfico muestra por qué es importante la variable de generación. Sin ella (panel de arriba), el modelo subestima los carros 2016 a 2018 y sobreestima los 2011 a 2015, con un salto justo en el cambio de generación. Con ella (panel de abajo), los residuos quedan alrededor de cero en todos los años.

par(mfrow = c(2, 2))
plot(m3_log)

par(mfrow = c(1, 1))

bp_log <- bptest(m3_log)
sw_log <- shapiro.test(residuals(m3_log))
dw_log <- dwtest(m3_log)
kable(tabla_pruebas(bp_log, sw_log, dw_log), digits = 3, caption = "Pruebas de supuestos para M3-log") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Pruebas de supuestos para M3-log
	Supuesto	Prueba	Estadistico	p_valor	Se_rechaza_al_5
BP	Varianza constante	Breusch-Pagan	10.422	0.0339	Sí
W	Normalidad de los errores	Shapiro-Wilk	0.848	<0.001	Sí
DW	Independencia	Durbin-Watson	1.904	0.1975	No

vif_log <- vif(m3_log)
vif_log

##        edad      km_mil transmision  generacion 
##    7.397329    4.351257    1.042872    4.064125

influencia <- data.frame(id = datos$id, modelo = datos$modelo, precio_m = datos$precio_m,
                         km_mil = datos$km_mil, transmision = datos$transmision,
                         residuo_est = rstudent(m3_log), apalancamiento = hatvalues(m3_log),
                         cook = cooks.distance(m3_log))

ggplot(influencia, aes(x = apalancamiento, y = residuo_est)) +
  geom_point(aes(size = cook), color = "steelblue", alpha = 0.6) +
  geom_hline(yintercept = c(-3, 3), linetype = "dashed") +
  geom_text(data = filter(influencia, abs(residuo_est) > 3), aes(label = id), size = 3.5, vjust = -1) +
  labs(title = "Residuos estudentizados frente a apalancamiento (M3-log)",
       x = "Apalancamiento", y = "Residuo estudentizado", size = "Distancia de Cook")

influencia %>%
  arrange(desc(cook)) %>%
  head(5) %>%
  kable(digits = 3, caption = "Anuncios más influyentes en M3-log") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Anuncios más influyentes en M3-log
	id	modelo	precio_m	km_mil	transmision	residuo_est	apalancamiento	cook
116	123	2017	28.9	12.345	Mecánica	-9.860	0.035	0.528
74	77	2011	22.0	102.000	Mecánica	-5.383	0.011	0.059
292	306	2010	26.0	176.000	Mecánica	-2.040	0.032	0.027
281	295	2012	29.0	145.000	Mecánica	-2.256	0.021	0.021
285	299	2009	25.8	129.000	Automática	-2.443	0.017	0.021

n_cook_log <- sum(influencia$cook > 4 / n_final)

Linealidad: la curva en U desapareció del gráfico de residuos contra valores ajustados: ahora los residuos quedan repartidos alrededor de cero.
Varianza constante: Breusch-Pagan da p = 0.034. Todavía queda algo de heterocedasticidad, aunque mucho menos que en el modelo lineal.
Normalidad: Shapiro-Wilk sigue rechazando (p < 0.001). El gráfico Q-Q muestra que el problema está en la cola izquierda: unos pocos anuncios tienen precios muy por debajo de lo esperado. Por eso en la sección 2.5 reviso si los intervalos de predicción contienen realmente el 95 % de los precios.
Multicolinealidad: el VIF más alto es 7.4 (edad). Es moderado y no afecta las predicciones.
Independencia: Durbin-Watson da 1.9.
Observaciones influyentes: 13 anuncios superan el umbral de 4/n. El más influyente es el id 123, el 2017 con “12.345 km” anunciado a 28,9 millones que ya había aparecido como atípico en el Punto 1, seguido del id 77, un 2011 anunciado a 22 millones. Los dos están muy por debajo de lo que el modelo espera para sus características. No hay pruebas de que sean errores, así que los mantengo; en la sección 2.6 muestro qué pasa si se excluyen. Un uso práctico del modelo sería justamente detectar anuncios como estos, que pueden ser errores o engaños.

MAE del modelo transformado frente al original

Para comparar con el modelo original, las predicciones de los modelos logarítmicos se devuelven a millones de pesos con la función exponencial, y el MAE se calcula en esa escala. Todos los modelos múltiples se evalúan con las mismas 100 particiones del Punto 1, y la selección por AIC se repite dentro de cada partición de entrenamiento:

f_m1 <- "precio_m ~ edad + km_mil + transmision"
f_m2 <- "precio_m ~ edad + km_mil + transmision + generacion + region + color_g"
f_m3 <- "precio_m ~ edad + km_mil + transmision + generacion"
a_log <- function(f) sub("precio_m", "log(precio_m)", f)   # misma fórmula con log(precio)

cv_res[["M1: lineal"]] <- cv_modelo(f_m1, particiones)
cv_res[["M2: lineal"]] <- cv_modelo(f_m2, particiones)
cv_res[["M2-AIC: lineal"]] <- cv_modelo(f_m2, particiones, paso_aic = TRUE)
cv_res[["M1-log"]] <- cv_modelo(a_log(f_m1), particiones, log = TRUE)
cv_res[["M2-log"]] <- cv_modelo(a_log(f_m2), particiones, log = TRUE)
cv_res[["M2-log-AIC"]] <- cv_modelo(a_log(f_m2), particiones, log = TRUE, paso_aic = TRUE)
cv_res[["M3-log"]] <- cv_modelo(a_log(f_m3), particiones, log = TRUE)

tabla_cv <- resumen_cv(cv_res)

tabla_cv %>%
  filter(Modelo %in% c("M2-AIC: lineal", "M2-log-AIC", "M3-log")) %>%
  kable(digits = 2, caption = "Modelo lineal frente a modelos transformados (errores en millones de pesos)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Modelo lineal frente a modelos transformados (errores en millones de pesos)
Modelo	MAE	DE_MAE	RMSE	MAPE	Cobertura_IP95	Mejora_vs_mediana
M2-log-AIC	2.87	0.51	4.24	5.84	96.02	80.50
M3-log	2.88	0.52	4.25	5.88	96.02	80.45
M2-AIC: lineal	3.50	0.57	4.91	7.46	94.95	76.18

mae_m2aic <- mae_de(tabla_cv, "M2-AIC: lineal")
mae_m2logaic <- mae_de(tabla_cv, "M2-log-AIC")
mae_m3 <- mae_de(tabla_cv, "M3-log")

Con la transformación, el MAE baja de 3.5 millones (M2-AIC lineal) a 2.87 millones (M2-log-AIC), una reducción de 18 %. M3-log, sin la región, tiene un MAE de 2.88 millones.

2.5 Comparación de todos los modelos y escenarios

kable(tabla_cv, digits = 2,
      caption = "Validación cruzada 10 x 10 con las mismas particiones (errores en millones de pesos; MAPE, cobertura y mejora en %)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE) %>%
  row_spec(which(tabla_cv$Modelo == "M3-log"), bold = TRUE, background = "#fde9d9")

Validación cruzada 10 x 10 con las mismas particiones (errores en millones de pesos; MAPE, cobertura y mejora en %)
Modelo	MAE	DE_MAE	RMSE	MAPE	Cobertura_IP95	Mejora_vs_mediana
M2-log	2.86	0.51	4.22	5.81	95.85	80.59
M2-log-AIC	2.87	0.51	4.24	5.84	96.02	80.50
M3-log	2.88	0.52	4.25	5.88	96.02	80.45
M1-log	3.48	0.51	4.62	7.13	97.29	76.38
M2-AIC: lineal	3.50	0.57	4.91	7.46	94.95	76.18
M2: lineal	3.56	0.56	4.94	7.58	95.55	75.77
M1: lineal	4.09	0.55	5.30	9.09	95.99	72.20
S1: edad	4.39	0.59	5.75	9.66	96.32	70.17
S2: kilometraje	7.21	1.23	9.24	16.18	94.98	51.01
S3: transmisión	14.40	1.48	16.71	31.32	96.89	2.13
Mediana	14.71	1.76	19.12	26.79	NA	0.00

tabla_cv %>%
  mutate(Modelo = factor(Modelo, levels = rev(Modelo)),
         Tipo = case_when(Modelo == "M3-log" ~ "Modelo elegido",
                          Modelo == "Mediana" ~ "Referencia",
                          TRUE ~ "Otros")) %>%
  ggplot(aes(x = MAE, y = Modelo, fill = Tipo)) +
  geom_col() +
  geom_text(aes(label = round(MAE, 2)), hjust = -0.2, size = 3.5) +
  scale_fill_manual(values = c("Modelo elegido" = "darkorange", "Referencia" = "gray60", "Otros" = "steelblue")) +
  scale_x_continuous(limits = c(0, max(tabla_cv$MAE) * 1.12)) +
  labs(title = "Error promedio de predicción de cada modelo",
       subtitle = "MAE en validación cruzada, millones de pesos",
       x = "MAE (millones de pesos)", y = NULL, fill = NULL) +
  theme(legend.position = "top")

Lo que muestra la comparación:

Los modelos múltiples superan a los simples. M3-log reduce el error de S1 de 4.39 a 2.88 millones y le gana en el 100 % de las particiones.
El logaritmo mejora la predicción en pesos. M3-log le gana al mejor modelo lineal (M2-AIC) en el 98 % de las particiones.
La generación es la variable que más aporta después de la edad. De M1-log a M3-log, que solo se diferencian en esa variable, el MAE baja de 3.48 a 2.88 millones.
La región y el color casi no mejoran la predicción. El modelo con menor MAE es M2-log (2.86 millones). M3-log está a 20 mil pesos de él, una diferencia muy pequeña frente a un error promedio de casi 3 millones. Cuando dos modelos predicen prácticamente igual, conviene quedarse con el más sencillo.
Los intervalos de predicción funcionan. El intervalo del 95 % de M3-log contiene el precio real en el 96 % de los casos de prueba, muy cerca del 95 % que promete.

Modelo elegido: M3-log, es decir, log(precio) en función de la edad, el kilometraje, la transmisión y la generación. Tiene un error promedio de 2.88 millones (5.9 % del precio), frente a 14.71 millones de la referencia. Además usa solo variables que cualquier vendedor conoce de su carro, sus coeficientes tienen una lectura sencilla en porcentajes y el patrón de depreciación tiene sentido.

Escenarios hipotéticos

Planteo dos casos realistas, dentro del rango de los datos:

escenarios <- data.frame(
  caso = c("1. Vendedora particular", "2. Compraventa (retoma)"),
  modelo = c(2015, 2019),
  km_mil = c(85, 45),
  transmision = factor(c("Mecánica", "Automática"), levels = levels(datos$transmision))
)
escenarios$edad <- 2022 - escenarios$modelo
escenarios$generacion <- factor(ifelse(escenarios$modelo >= 2016, "Nueva (2016 en adelante)", "Anterior (2008 a 2015)"),
                                levels = levels(datos$generacion))

pred_esc <- exp(predict(m3_log, newdata = escenarios, interval = "prediction", level = 0.95))

comparables_med <- sapply(1:2, function(i) median(datos$precio_m[datos$modelo == escenarios$modelo[i] &
                                                                  datos$transmision == escenarios$transmision[i]]))
comparables_n <- sapply(1:2, function(i) sum(datos$modelo == escenarios$modelo[i] &
                                             datos$transmision == escenarios$transmision[i]))

data.frame(
  Caso = escenarios$caso,
  Modelo = escenarios$modelo,
  Km = escenarios$km_mil * 1000,
  Transmision = escenarios$transmision,
  Prediccion = pred_esc[, "fit"],
  IP95_inf = pred_esc[, "lwr"],
  IP95_sup = pred_esc[, "upr"],
  Mediana_comparables = comparables_med,
  N_comparables = comparables_n
) %>%
  kable(digits = 1, caption = "Predicciones de M3-log en millones de pesos con intervalo de predicción del 95 %") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Predicciones de M3-log en millones de pesos con intervalo de predicción del 95 %
Caso	Modelo	Km	Transmision	Prediccion	IP95_inf	IP95_sup	Mediana_comparables	N_comparables
Vendedora particular	2015	85000	Mecánica	42.5	36.0	50.1	43.9	11
2. Compraventa (retoma)	2019	45000	Automática	68.0	57.7	80.1	70.0	12

# contraste: mismo carro, un año de diferencia entre generaciones
contraste <- data.frame(modelo = c(2015, 2016), km_mil = 70,
                        transmision = factor("Automática", levels = levels(datos$transmision)))
contraste$edad <- 2022 - contraste$modelo
contraste$generacion <- factor(ifelse(contraste$modelo >= 2016, "Nueva (2016 en adelante)", "Anterior (2008 a 2015)"),
                               levels = levels(datos$generacion))
pred_contraste <- exp(predict(m3_log, newdata = contraste))
pred_contraste

##        1        2 
## 45.44114 56.72142

Caso 1. Una vendedora particular quiere publicar su Mazda2 2015 mecánico con 85.000 km. El modelo estima un precio de publicación de 42.5 millones, con un intervalo del 95 % entre 36 y 50.1 millones. En la base hay 11 anuncios de 2015 mecánicos, con una mediana de 43.9 millones, así que la predicción es coherente con el mercado. Si ella quiere vender rápido, le conviene publicar cerca del centro del intervalo. Si el carro está en muy buen estado y tiene cómo demostrarlo, puede ubicarse en la parte alta. Por encima de 50.1 millones estaría pidiendo más que lo que piden casi todos los vendedores de carros parecidos.

Caso 2. Una compraventa recibe en parte de pago un Mazda2 2019 automático con 45.000 km. El modelo estima un precio anunciado de 68 millones, con un intervalo entre 57.7 y 80.1 millones (los 12 anuncios de 2019 automáticos tienen una mediana de 70 millones). Como el modelo estima precios de anuncio y no de venta, para la compraventa esa cifra es más bien un techo de reventa: a eso tiene que descontarle su margen, el alistamiento y lo que encuentre en el peritaje. Una oferta prudente estaría por debajo del límite inferior del intervalo.

Contraste entre generaciones. Un Mazda2 automático 2015 con 70.000 km se estima en 45.4 millones, y el mismo carro modelo 2016 en 56.7 millones. Un solo año de diferencia representa 11.3 millones, porque el 2016 ya es de la generación nueva. Un modelo sin esa variable no podría ver esa diferencia.

2.6 Verificaciones adicionales

En esta sección reviso si las conclusiones cambian con otra división al azar de los datos o con otras decisiones de limpieza, y si los intervalos de predicción funcionan.

Otras semillas

sens_multiples <- data.frame()
for (s in semillas) {
  part_s <- crear_particiones(s, R = 5)
  sens_multiples <- rbind(sens_multiples, data.frame(
    Semilla = s,
    Mediana = sens_simples$Mediana[sens_simples$Semilla == s],
    S1 = sens_simples$S1[sens_simples$Semilla == s],
    M2_AIC_lineal = mean(mae_por_particion(cv_modelo(f_m2, part_s, paso_aic = TRUE))),
    M3_log = mean(mae_por_particion(cv_modelo(a_log(f_m3), part_s, log = TRUE)))
  ))
}
kable(sens_multiples, digits = 3, caption = "MAE con cinco semillas distintas (millones de pesos)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

MAE con cinco semillas distintas (millones de pesos)
Semilla	Mediana	S1	M2_AIC_lineal	M3_log
1	14.712	4.385	3.488	2.863
42	14.711	4.392	3.489	2.871
123	14.713	4.393	3.514	2.869
777	14.698	4.394	3.479	2.875
2024	14.706	4.396	3.495	2.864

orden_multiples <- all(sens_multiples$M3_log < sens_multiples$M2_AIC_lineal &
                       sens_multiples$M2_AIC_lineal < sens_multiples$S1 &
                       sens_multiples$S1 < sens_multiples$Mediana)

Con las cinco semillas, el MAE de M3-log va de 2.86 a 2.88 millones. El orden entre modelos es siempre el mismo: M3-log, M2-AIC lineal, S1 y la mediana.

Decisiones de limpieza

crudos <- mazda %>%
  mutate(precio_m = precio / 1e6,
         edad = 2022 - modelo,
         km_mil = kilometraje / 1000,
         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)")))

bases <- list(
  "Base limpia (principal)" = datos,
  "Excluyendo los km corregidos" = datos %>% filter(!km_en_miles),
  "Excluyendo residuos estudentizados > 3" = datos[abs(influencia$residuo_est) <= 3, ],
  "Base original sin limpiar" = crudos
)

limpieza <- data.frame()
for (nombre in names(bases)) {
  m <- lm(log(precio_m) ~ edad + km_mil + transmision + generacion, data = bases[[nombre]])
  limpieza <- rbind(limpieza, data.frame(
    Base = nombre, n = nrow(bases[[nombre]]),
    Edad_pct = ef(coef(m)[2]), Km10_pct = ef(10 * coef(m)[3]),
    Automatica_pct = ef(coef(m)[4]), Generacion_pct = ef(coef(m)[5]),
    Error_residual = sigma(m), R2 = summary(m)$r.squared, row.names = NULL))
}
kable(limpieza, digits = 3, caption = "M3-log estimado con distintas versiones de la base") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

M3-log estimado con distintas versiones de la base
Base	n	Edad_pct	Km10_pct	Automatica_pct	Generacion_pct	Error_residual	R2
Base limpia (principal)	299	-5.468	-0.488	6.108	17.999	0.083	0.940
Excluyendo los km corregidos	288	-5.438	-0.559	6.285	17.870	0.084	0.939
Excluyendo residuos estudentizados > 3	296	-5.078	-0.794	4.936	19.486	0.066	0.961
Base original sin limpiar	314	-4.676	-0.166	7.793	26.701	0.193	0.747

veces_error <- limpieza$Error_residual[4] / limpieza$Error_residual[1]

Corregir o excluir los kilometrajes escritos en miles da casi los mismos resultados, así que esa decisión no cambia nada importante. Excluir los atípicos extremos reduce el error del modelo, pero no cambia las conclusiones. En cambio, no limpiar la base sí importa: con los 314 anuncios originales el error residual es 2.3 veces mayor y el efecto de la generación se infla, porque las camionetas CX-5 y CX-9 de 160 a 170 millones quedan clasificadas como Mazda2 de la generación nueva.

Calibración y error por segmentos

oof <- cv_res[["M3-log"]] %>%
  left_join(datos %>% select(id, generacion, transmision, region, edad), by = "id") %>%
  mutate(rango_edad = cut(edad, breaks = c(-Inf, 2, 6, 10, Inf),
                          labels = c("0 a 2 años", "3 a 6 años", "7 a 10 años", "11 a 14 años")))

por_segmento <- function(variable) {
  oof %>%
    group_by(Grupo = as.character(.data[[variable]])) %>%
    summarise(Predicciones = n(),
              MAE = mean(abs(real - pred)),
              MAPE = 100 * mean(abs(real - pred) / real),
              Cobertura_IP95 = 100 * mean(real >= inf & real <= sup),
              Sesgo = mean(pred - real)) %>%
    mutate(Variable = variable) %>%
    select(Variable, everything())
}
segmentos <- bind_rows(por_segmento("generacion"), por_segmento("transmision"),
                       por_segmento("rango_edad"), por_segmento("region"))
kable(segmentos, digits = 2, caption = "Desempeño de M3-log por segmento (predicciones fuera de muestra)") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Desempeño de M3-log por segmento (predicciones fuera de muestra)
Variable	Grupo	Predicciones	MAE	MAPE	Cobertura_IP95	Sesgo
generacion	Anterior (2008 a 2015)	1670	1.89	5.46	97.01	-0.06
generacion	Nueva (2016 en adelante)	1320	4.12	6.40	94.77	-0.35
transmision	Automática	1640	3.03	5.23	96.46	-0.43
transmision	Mecánica	1350	2.69	6.66	95.48	0.10
rango_edad	0 a 2 años	500	5.23	6.62	94.60	-0.89
rango_edad	11 a 14 años	730	2.01	6.48	95.89	-0.21
rango_edad	3 a 6 años	820	3.45	6.27	94.88	-0.02
rango_edad	7 a 10 años	940	1.80	4.67	97.87	0.06
region	Antioquia	800	3.32	6.24	92.88	-1.33
region	Bogotá y Cundinamarca	550	3.25	7.29	98.18	0.95
region	Eje Cafetero	430	2.27	4.67	100.00	-0.93
region	Otros departamentos	710	2.22	4.84	98.59	0.34
region	Valle del Cauca	500	3.22	6.24	91.60	0.28

sesgo_ant <- segmentos$Sesgo[segmentos$Grupo == "Antioquia"]
sesgo_bog <- segmentos$Sesgo[segmentos$Grupo == "Bogotá y Cundinamarca"]

oof %>%
  group_by(id, generacion) %>%
  summarise(real = first(real), pred = mean(pred), inf = mean(inf), sup = mean(sup), .groups = "drop") %>%
  ggplot(aes(x = pred, y = real, color = generacion)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_errorbar(aes(xmin = inf, xmax = sup), width = 0, alpha = 0.2, orientation = "y") +
  geom_point(alpha = 0.7) +
  scale_color_manual(values = c("steelblue", "darkorange")) +
  labs(title = "Precio real frente a predicción fuera de muestra (M3-log)",
       subtitle = "Barras horizontales: intervalo de predicción del 95 %",
       x = "Predicción (millones de pesos)", y = "Precio anunciado (millones de pesos)", color = "Generación") +
  theme(legend.position = "top")

El error relativo es parecido entre segmentos (el MAPE va de 4.7 % a 7.3 %) y la cobertura de los intervalos va de 91.6 % a 100 %. En pesos, el error es mayor en los carros nuevos, porque son más caros y porque en ellos pesa más la versión, que el modelo no conoce. Por región aparece un sesgo moderado: el modelo tiende a subestimar los precios de Antioquia (-1.33 millones en promedio) y a sobreestimar los de Bogotá y Cundinamarca (0.95 millones). Es coherente con lo que encontró el AIC, pero es pequeño frente al error total. Si el modelo se fuera a usar en una sola ciudad, convendría recalibrarlo con datos de esa ciudad.

Resumen de las verificaciones

data.frame(
  Verificacion = c("Cinco semillas de validación cruzada", "Versiones de la limpieza",
                   "Cobertura del intervalo del 95 %"),
  Resultado = c(
    paste0("MAE de M3-log entre ", round(min(sens_multiples$M3_log), 2), " y ", round(max(sens_multiples$M3_log), 2),
           " millones; orden de los modelos ", ifelse(orden_multiples, "igual en todas", "casi igual")),
    paste0("Corregir los km o quitar los atípicos no cambia las conclusiones; sin limpieza el error se multiplica por ",
           round(veces_error, 1)),
    paste0(round(cob_m3, 1), " % fuera de muestra; entre ", round(min(segmentos$Cobertura_IP95)), " % y ",
           round(max(segmentos$Cobertura_IP95)), " % por segmento")
  )
) %>%
  kable(caption = "Resumen de las verificaciones") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Resumen de las verificaciones
Verificacion	Resultado
Cinco semillas de validación cruzada	MAE de M3-log entre 2.86 y 2.88 millones; orden de los modelos igual en todas
Versiones de la limpieza	Corregir los km o quitar los atípicos no cambia las conclusiones; sin limpieza el error se multiplica por 2.3
Cobertura del intervalo del 95 %	96 % fuera de muestra; entre 92 % y 100 % por segmento

Punto 3. Análisis crítico

¿Quién utilizaría el modelo y para qué decisión?

Personas que venden un Mazda2 usado, para decidir el precio de publicación. Es el uso más directo, porque el modelo se entrenó justamente con precios de publicación.
Compradores, para saber si un anuncio está caro o barato frente a carros parecidos y tener un argumento para negociar.
Compraventas y concesionarios, como primera valoración de un carro que reciben en parte de pago, antes de pagar un peritaje.
La plataforma de anuncios, para sugerir un rango de precio al publicar y marcar anuncios con precios muy por fuera de lo esperado, como el Mazda2 2022 a 6,8 millones o el 2017 con “12.345 km” a mitad de precio, que pueden ser errores o engaños.

¿Qué significa el error de predicción en la práctica? ¿Sería aceptable?

En promedio, el modelo se equivoca en 2.88 millones, que equivale al 5.9 % del precio. La mitad de los errores son menores a 1.98 millones. El intervalo del 95 % tiene un ancho promedio de más o menos 17 %: para un carro de 60 millones, el rango iría aproximadamente de 50 a 70 millones.

Para fijar un precio de publicación o para decidir si vale la pena mirar un anuncio, este error me parece aceptable: es 5.1 veces menor que el de la referencia y los intervalos están bien calibrados. Como valor definitivo para una compra no es suficiente, porque una diferencia de 3 a 5 millones es mucha plata para cualquiera de las dos partes. Ahí el modelo debe acompañarse de una revisión del carro.

¿Qué consecuencias tendría sobreestimar o subestimar?

Para quien vende: si el modelo sobreestima, publica muy caro, el carro no se vende y le toca bajar el precio después de perder tiempo. Si subestima, vende rápido pero pierde plata.
Para una compraventa: sobreestimar es lo más grave, porque pagaría de más por el carro y perdería su margen al revenderlo. Subestimar le haría perder el negocio frente a otra compraventa. Como el costo no es el mismo en ambos casos, le conviene usar el límite inferior del intervalo como referencia para comprar.
Para la plataforma: un umbral muy estricto marcaría como sospechosos anuncios normales y molestaría a los usuarios; uno muy flexible dejaría pasar anuncios engañosos.

En la validación cruzada el modelo sobreestima en el 47 % de los casos, así que no tiene un sesgo sistemático hacia arriba ni hacia abajo.

¿Cuáles son sus principales limitaciones y cuándo no debería utilizarse?

Precio anunciado no es precio de venta. El modelo reproduce lo que piden los vendedores. Lo más probable es que el precio final de venta sea menor.
Los datos son de 2022. Desde entonces los precios cambiaron con la inflación y con el mercado de carros nuevos. Antes de usar el modelo hoy habría que volver a estimarlo con anuncios recientes. Los porcentajes de depreciación probablemente se mantienen mejor que los valores en pesos.
Faltan variables importantes: la versión, el estado del carro, los accidentes, el número de dueños, si está blindado y el historial de mantenimiento. Buena parte del error que queda se debe a ellas.
La muestra es pequeña y de una sola plataforma: unos 300 anuncios de OLX en unas pocas semanas. Algunos concesionarios publican varios carros, lo que hace que las observaciones no sean del todo independientes.
No debería usarse para carros fuera del rango de los datos (antes de 2008, después de 2022 o con más de 200.000 km), para carros de salvamento, chocados, blindados o modificados, para otras referencias o marcas, ni como avalúo oficial (para seguros y créditos existe la guía de valores de Fasecolda). Tampoco sirve para sacar conclusiones causales.

¿Recomiendo aplicarlo? ¿Qué información adicional necesitaría?

Sí recomiendo aplicarlo, pero como herramienta de referencia y de alerta, no como precio final, y con tres condiciones:

Volver a estimarlo con anuncios actuales antes de usarlo.
Mostrar siempre el rango de precio y no solo un número.
En decisiones de compra, complementarlo con un peritaje.

Para mejorarlo necesitaría precios de venta reales y no solo de anuncio, la versión del carro (Prime, Touring o Grand Touring), que hoy no está en la base, la fecha de publicación y los días que duró el anuncio, información sobre el estado del carro y su historial, y anuncios de otras plataformas para tener una muestra más representativa.

Conclusiones

El modelo elegido estima el logaritmo del precio con la edad, el kilometraje, la transmisión y la generación. Su error promedio fuera de muestra es de 2.88 millones (5.9 % del precio), frente a 14.71 millones si se predijera siempre la mediana, y sus intervalos del 95 % contienen el precio real en el 96 % de los casos.
Cada año de antigüedad resta en promedio 5.5 % del precio, la caja automática suma 6.1 % y la generación nueva suma 18 %. Una vez se conoce la edad, el kilometraje aporta poco.
La transformación logarítmica corrigió la curvatura y redujo el error frente al mejor modelo lineal (de 3.5 a 2.88 millones).
Los resultados se mantienen con otras semillas y sin los anuncios más extremos. La limpieza de la base fue la decisión que más impacto tuvo.
El modelo es útil como referencia de precio de publicación y para detectar anuncios con precios anómalos, siempre que se actualice con datos recientes, se comunique como un rango y no se use como precio de venta definitivo.

Referencias

Autocosmos (27 de marzo de 2015). Mazda2 2016 llega a Colombia desde 43 millones de pesos. https://noticias.autocosmos.com.co/2015/03/27/mazda2-2016-llega-a-colombia-desde-43-millones-de-pesos

Información de la sesión

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.1 knitr_1.52       ggplot2_4.0.3    tidyr_1.3.2     
##  [5] dplyr_1.2.1      lmtest_0.9-40    zoo_1.9-0        car_3.1-5       
##  [9] carData_3.0-6    MASS_7.3-66      readxl_1.6.0    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     xml2_1.6.0         stringi_1.8.9     
##  [5] lattice_0.22-9     digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.1         RColorBrewer_1.1-3 fastmap_1.2.0      Matrix_1.7-5      
## [13] cellranger_1.1.0   jsonlite_2.0.0     Formula_1.2-6      mgcv_1.9-4        
## [17] purrr_1.2.2        viridisLite_0.4.3  scales_1.4.0       textshaping_1.0.5 
## [21] jquerylib_0.1.4    abind_1.4-8        cli_3.6.6          rlang_1.3.0       
## [25] splines_4.6.1      withr_3.0.3        cachem_1.1.0       yaml_2.3.12       
## [29] otel_0.2.0         tools_4.6.1        vctrs_0.7.3        R6_2.6.1          
## [33] stats4_4.6.1       lifecycle_1.0.5    stringr_1.6.0      pkgconfig_2.0.3   
## [37] pillar_1.11.1      bslib_0.12.0       gtable_0.3.6       glue_1.8.1        
## [41] systemfonts_1.3.2  xfun_0.60          tibble_3.3.1       tidyselect_1.2.1  
## [45] rstudioapi_0.19.0  farver_2.1.2       nlme_3.1-169       htmltools_0.5.9   
## [49] labeling_0.4.3     rmarkdown_2.32     svglite_2.2.2      compiler_4.6.1    
## [53] S7_0.2.2

Caso 2: Modelo de regresión lineal simple y múltiple

¿Cuánto pedir por un Mazda2 usado? Base 1: anuncios de OLX

Juan Esteban Vidal

11 de noviembre de 2026