Taller 2 - Modelos Lineales
California Housing
Juan Sánchez
12 de junio de 2026
1 Preparación
El archivo de Kaggle trae nombres distintos a los del taller, así que lo primero es dejar las cinco variables que vamos a estudiar con un significado claro: el valor medio de la vivienda, el ingreso medio del distrito, la antigüedad de las viviendas, y dos medidas por hogar que hay que construir —las habitaciones promedio y los ocupantes promedio—. El valor de la vivienda lo dejo en dólares y el ingreso en decenas de miles de dólares, tal como vienen en el archivo: no transformo nada, para mantener la coherencia con la naturaleza original de los datos.
# Manejo de datos y gráficos
library(tidyverse)
# Gráficos de correlación
library(corrplot)
library(PerformanceAnalytics)
# Asimetría y curtosis (coeficientes de forma para evaluar normalidad)
library(moments)
# Pruebas de supuestos: Breusch-Pagan (homocedasticidad) y Durbin-Watson (independencia)
library(lmtest)
# VIF para medir multicolinealidad
library(car)
# Pruebas de normalidad para muestras grandes (Lilliefors)
library(nortest)
# Formato de ejes y cifras en dólares
library(scales)
library(gridExtra)
# Salidas ordenadas de los modelos
library(broom)crudo <- read_csv("housing.csv")
# Me quedo solo con las cinco variables del taller y construyo las dos que no vienen
# directas (las medidas "por hogar" se obtienen dividiendo entre el número de hogares).
datos <- crudo %>%
transmute(
MedHouseVal = median_house_value, # valor medio de la vivienda (USD)
MedInc = median_income, # ingreso medio (decenas de miles USD)
HouseAge = housing_median_age, # antigüedad media de la vivienda (años)
AveRooms = total_rooms / households, # habitaciones promedio por hogar
AveOccup = population / households # ocupantes promedio por hogar
)
glimpse(datos)## Rows: 20,640
## Columns: 5
## $ MedHouseVal <dbl> 452600, 358500, 352100, 341300, 342200, 269700, 299200, 24…
## $ MedInc <dbl> 8.3252, 8.3014, 7.2574, 5.6431, 3.8462, 4.0368, 3.6591, 3.…
## $ HouseAge <dbl> 41, 21, 52, 52, 52, 52, 52, 52, 42, 52, 52, 52, 52, 52, 52…
## $ AveRooms <dbl> 6.984127, 6.238137, 8.288136, 5.817352, 6.281853, 4.761658…
## $ AveOccup <dbl> 2.555556, 2.109842, 2.802260, 2.547945, 2.181467, 2.139896…2 Análisis exploratorio
El objetivo natural de este conjunto de datos es explicar y, en lo posible, predecir el valor de las viviendas de un distrito a partir de sus características socioeconómicas, sobre todo del ingreso de sus habitantes. Esa es la pregunta que guía todo el análisis.
# Tamaño del problema: cuántos distritos (filas) y cuántas variables trae el archivo.
n_obs <- nrow(datos)
n_vars <- ncol(crudo) # variables del dataset original, antes de seleccionar
n_obs## [1] 20640n_vars## [1] 10El conjunto original reúne 20640 distritos (las observaciones) descritos por 10 variables. De todas ellas nos quedamos con las cinco que pide el taller.
# Resumen numérico de cada variable. Incluyo el centro (media y mediana), la dispersión
# absoluta (sd) y, sobre todo, el coeficiente de variación (cv = sd/media): es la
# dispersión RELATIVA y es lo único que permite comparar variables de escalas distintas.
resumen <- datos %>%
pivot_longer(everything(), names_to = "variable", values_to = "valor") %>%
group_by(variable) %>%
summarise(
media = mean(valor),
mediana = median(valor),
sd = sd(valor),
cv = sd(valor) / mean(valor),
min = min(valor),
max = max(valor),
.groups = "drop"
)
knitr::kable(resumen, digits = 2, caption = "Estadísticos descriptivos")| variable | media | mediana | sd | cv | min | max |
|---|---|---|---|---|---|---|
| AveOccup | 3.07 | 2.82 | 10.39 | 3.38 | 0.69 | 1243.33 |
| AveRooms | 5.43 | 5.23 | 2.47 | 0.46 | 0.85 | 141.91 |
| HouseAge | 28.64 | 29.00 | 12.59 | 0.44 | 1.00 | 52.00 |
| MedHouseVal | 206855.82 | 179700.00 | 115395.62 | 0.56 | 14999.00 | 500001.00 |
| MedInc | 3.87 | 3.53 | 1.90 | 0.49 | 0.50 | 15.00 |
El valor medio de las viviendas ronda los $206,856, con una mediana algo menor ($179,700); esa diferencia entre media y mediana ya anticipa una distribución estirada hacia los precios altos.
Sobre la dispersión hay que tener cuidado. Si nos guiamos por la desviación estándar en bruto, parecería que el valor de la vivienda es la variable más dispersa, pero eso es solo un efecto de su escala: está en dólares, mientras las demás se miden en años o en conteos por hogar. Comparando de forma justa con el coeficiente de variación, la variable que más varía en términos relativos es la ocupación por hogar, es decir, cuánta gente vive en promedio en cada vivienda.
# Distribución de cada variable EN SU ESCALA REAL. Uso facetas con escala libre en vez de
# estandarizar (z-score): así cada caja se ve en sus propias unidades —más honesto
# visualmente— y de paso quedan a la vista los valores atípicos de cada variable.
etiquetas <- c(
MedHouseVal = "Valor de vivienda (USD)",
MedInc = "Ingreso medio (decenas de miles USD)",
HouseAge = "Antigüedad de la vivienda (años)",
AveRooms = "Habitaciones por hogar",
AveOccup = "Ocupantes por hogar"
)
datos %>%
pivot_longer(everything(), names_to = "variable", values_to = "valor") %>%
ggplot(aes(y = valor)) +
geom_boxplot(fill = "steelblue", alpha = 0.6, outlier.size = 0.4) +
facet_wrap(~ variable, scales = "free", nrow = 2,
labeller = as_labeller(etiquetas)) +
labs(title = "Distribución de cada variable en su escala original",
y = NULL) +
theme_minimal() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
Valores atípicos hay, y bastantes. Los diagramas de caja muestran colas largas en las habitaciones por hogar y, sobre todo, en la ocupación por hogar, donde aparecen distritos con cifras extremas. A eso se suma un detalle que conviene no perder de vista: los datos están topados. Hay 965 viviendas pegadas al límite de $500,001 y 51 distritos en el ingreso máximo de 15. No elimino estas observaciones —prefiero trabajar con el conjunto tal como viene—, pero las tendré presentes al interpretar los modelos, porque son justamente las que más tarde tensionan los supuestos.
Para ver mejor dónde están esos extremos, grafico el precio frente a las tres variables más relacionadas con él. Cada punto es un distrito; dejo el eje horizontal libre en cada panel para que se note la escala real de cada relación.
# Dispersión del precio frente a sus variables más asociadas. El eje x libre deja ver cómo
# unos pocos distritos con ocupación o habitaciones por hogar disparatadas estiran la
# escala y separan la nube principal de esos puntos extremos.
datos %>%
select(MedHouseVal, MedInc, AveRooms, AveOccup) %>%
pivot_longer(-MedHouseVal, names_to = "variable", values_to = "valor") %>%
ggplot(aes(valor, MedHouseVal)) +
geom_point(alpha = 0.15, size = 0.5, color = "steelblue") +
facet_wrap(~ variable, scales = "free_x",
labeller = as_labeller(etiquetas)) +
scale_y_continuous(labels = dollar) +
labs(title = "Precio vs. variables más relacionadas",
x = NULL, y = "Valor de vivienda (USD)") +
theme_minimal()
El gráfico del ingreso muestra la nube ordenada y creciente que ya esperábamos. Los de habitaciones y ocupación por hogar, en cambio, aparecen casi pegados al eje vertical: toda la masa de datos se concentra a la izquierda y un puñado de distritos con valores absurdamente altos estira la escala hacia la derecha. Conviene medir cuántos son.
# Cuento los valores extremos con la regla del rango intercuartílico (fuera de
# Q1 - 1.5*IQR o Q3 + 1.5*IQR) y los comparo con el total de distritos, para dimensionar
# el problema frente al volumen de datos.
n <- nrow(datos)
extremos <- function(x) {
q <- quantile(x, c(.25, .75)); i <- diff(q)
sum(x < q[1] - 1.5 * i | x > q[2] + 1.5 * i)
}
ext_occup <- extremos(datos$AveOccup)
ext_rooms <- extremos(datos$AveRooms)
tibble(
variable = c("Ocupantes por hogar", "Habitaciones por hogar"),
extremos = c(ext_occup, ext_rooms),
porcentaje = percent(c(ext_occup, ext_rooms) / n, accuracy = 0.1),
maximo = round(c(max(datos$AveOccup), max(datos$AveRooms)), 1)
)| variable | extremos | porcentaje | maximo |
|---|---|---|---|
| Ocupantes por hogar | 711 | 3.4% | 1243.3 |
| Habitaciones por hogar | 511 | 2.5% | 141.9 |
En volumen son pocos —alrededor del 3.4% de los distritos en la ocupación por hogar y el 2.5% en las habitaciones—, pero sus valores no tienen sentido físico: hasta 1243 personas o 142 habitaciones por hogar son cifras imposibles, casi seguro errores de registro o distritos atípicos (zonas con muy pocos hogares censados). Aunque sean una fracción mínima del conjunto, pesan mucho en las medidas que se basan en la media, como la correlación de Pearson. Antes de cerrar un modelo definitivo, lo recomendable sería darles un tratamiento coherente con el contexto —por ejemplo recortar (winsorizar) o filtrar los valores físicamente imposibles, o usar métodos robustos— y volver a evaluar el modelo ya sin esa distorsión; es muy probable que así mejore tanto el ajuste como el cumplimiento de los supuestos.
3 Covarianzas y correlaciones
# Matriz de covarianzas: mide si dos variables se mueven juntas, pero en unidades crudas,
# así que su tamaño depende de la escala de cada variable (no es comparable entre pares).
cov_mat <- cov(datos)
knitr::kable(cov_mat, digits = 2, caption = "Matriz de covarianzas")| MedHouseVal | MedInc | HouseAge | AveRooms | AveOccup | |
|---|---|---|---|---|---|
| MedHouseVal | 1.331615e+10 | 150847.48 | 153398.80 | 43382.56 | -28449.40 |
| MedInc | 1.508475e+05 | 3.61 | -2.85 | 1.54 | 0.37 |
| HouseAge | 1.533988e+05 | -2.85 | 158.40 | -4.77 | 1.72 |
| AveRooms | 4.338256e+04 | 1.54 | -4.77 | 6.12 | -0.12 |
| AveOccup | -2.844940e+04 | 0.37 | 1.72 | -0.12 | 107.87 |
La mayor covarianza positiva aparece entre el valor de la vivienda y el ingreso medio. Conviene leerla con cautela: como la covarianza depende de la escala, cualquier par que incluya el valor de la vivienda (en dólares) va a dominar la matriz por su tamaño y no necesariamente por la fuerza de la relación. Por eso, para concluir bien, es mejor pasar a las correlaciones. También hay relaciones negativas: la ocupación por hogar se mueve en sentido contrario al valor de la vivienda y al ingreso —donde viven más personas por hogar, las viviendas valen menos y los ingresos son más bajos—.
Antes de calcular las correlaciones conviene decidir qué coeficiente usar, y eso depende de si las variables son normales. La de Pearson mide relación lineal pero supone normalidad y es muy sensible a los valores atípicos; la de Spearman trabaja con rangos, así que capta relaciones monótonas y aguanta bien la no-normalidad y los extremos. Por eso primero reviso la normalidad de cada variable.
# Pruebo normalidad teniendo en cuenta el volumen: Shapiro está limitado a 5000 datos, así
# que uso Anderson-Darling (admite n grande). Como con ~20.000 observaciones casi cualquier
# prueba rechaza, acompaño el p-valor con la asimetría, que mide la magnitud real del sesgo.
normalidad <- datos %>%
pivot_longer(everything(), names_to = "variable", values_to = "valor") %>%
group_by(variable) %>%
summarise(asimetria = skewness(valor),
ad_pvalor = ad.test(valor)$p.value,
.groups = "drop")
knitr::kable(normalidad, digits = c(0, 2, 4),
caption = "Normalidad de cada variable (Anderson-Darling)")| variable | asimetria | ad_pvalor |
|---|---|---|
| AveOccup | 97.63 | 0 |
| AveRooms | 20.70 | 0 |
| HouseAge | 0.06 | 0 |
| MedHouseVal | 0.98 | 0 |
| MedInc | 1.65 | 0 |
Ninguna variable pasa la prueba (todos los p-valores son prácticamente cero) y algunas están fuertemente sesgadas: la asimetría de las habitaciones por hogar y, sobre todo, de la ocupación por hogar es enorme, justo por esos distritos atípicos que vimos en los boxplots. Con este panorama, quedarse solo con Pearson sería arriesgado, así que reporto los dos coeficientes y los comparo.
# La correlación es la covarianza "normalizada" (entre -1 y 1), ya comparable entre pares.
# Calculo Pearson (relación lineal) y Spearman (rangos, robusto a la no-normalidad).
# Para la multicolinealidad uso solo las correlaciones ENTRE predictores (sin el precio).
cor_mat <- cor(datos) # Pearson
cor_spearman <- cor(datos, method = "spearman") # Spearman
knitr::kable(cor_mat, digits = 3, caption = "Correlación de Pearson (lineal)")| MedHouseVal | MedInc | HouseAge | AveRooms | AveOccup | |
|---|---|---|---|---|---|
| MedHouseVal | 1.000 | 0.688 | 0.106 | 0.152 | -0.024 |
| MedInc | 0.688 | 1.000 | -0.119 | 0.327 | 0.019 |
| HouseAge | 0.106 | -0.119 | 1.000 | -0.153 | 0.013 |
| AveRooms | 0.152 | 0.327 | -0.153 | 1.000 | -0.005 |
| AveOccup | -0.024 | 0.019 | 0.013 | -0.005 | 1.000 |
knitr::kable(cor_spearman, digits = 3, caption = "Correlación de Spearman (rangos)")| MedHouseVal | MedInc | HouseAge | AveRooms | AveOccup | |
|---|---|---|---|---|---|
| MedHouseVal | 1.000 | 0.677 | 0.075 | 0.263 | -0.257 |
| MedInc | 0.677 | 1.000 | -0.147 | 0.644 | -0.044 |
| HouseAge | 0.075 | -0.147 | 1.000 | -0.231 | -0.025 |
| AveRooms | 0.263 | 0.644 | -0.231 | 1.000 | 0.019 |
| AveOccup | -0.257 | -0.044 | -0.025 | 0.019 | 1.000 |
cor_pred <- cor(select(datos, -MedHouseVal))
cor_pred_max <- max(abs(cor_pred[lower.tri(cor_pred)]))
corrplot(cor_mat, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", tl.srt = 45,
number.cex = 0.8)
# Vista combinada: en la diagonal el histograma de cada variable, abajo los diagramas de
# dispersión y arriba el valor de la correlación. Sirve para ver forma y relación a la vez.
chart.Correlation(datos, histogram = TRUE, pch = 19)
La variable más correlacionada con el valor de la vivienda es el ingreso medio, y en esto los dos coeficientes coinciden: Pearson da 0.688 y Spearman 0.677, una relación positiva y de fuerza moderada-alta —a mayor ingreso, viviendas más caras—. Que ambos coincidan es buena señal: la relación es sólida y no la sostiene un puñado de datos extremos.
Donde sí difieren es revelador. Para la ocupación por hogar, Pearson da casi cero (-0.024), como si no tuviera relación con el precio, pero Spearman sube a -0.257: los distritos con valores extremos de ocupación aplastan la correlación lineal, y solo al trabajar con rangos aparece la relación negativa real (más hacinamiento, viviendas más baratas). Algo parecido pasa con las habitaciones por hogar (0.152 en Pearson frente a 0.263 en Spearman). Es justo el aviso que adelantábamos: con datos no normales y con atípicos, Pearson se queda corto y Spearman da una lectura más fiel.
Entre los predictores, en cambio, las correlaciones son bajas —la mayor en valor absoluto es apenas 0.327—, de modo que no hay indicios de que se solapen entre sí. Esto importa porque una correlación alta entre predictores (multicolinealidad) volvería los coeficientes inestables, inflaría sus errores estándar y haría casi imposible separar el efecto propio de cada variable; aquí, por suerte, no es el caso.
4 Regresión lineal simple
Como el ingreso medio es la variable más asociada al precio, el modelo simple más razonable usa el ingreso para explicar el valor de la vivienda. Podrían plantearse otras regresiones simples (una por cada predictor), pero serían peores: el resto de variables tiene correlaciones mucho más débiles con el precio, así que el ingreso es la mejor elección para un único predictor.
# Ajuste por mínimos cuadrados de valor ~ ingreso. summary() entrega los coeficientes,
# sus errores estándar, los valores-p y el R2, que es lo que iré interpretando.
modelo_simple <- lm(MedHouseVal ~ MedInc, data = datos)
summary(modelo_simple)##
## Call:
## lm(formula = MedHouseVal ~ MedInc, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -540697 -55950 -16979 36978 434023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45085.6 1322.9 34.08 <2e-16 ***
## MedInc 41793.8 306.8 136.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 83740 on 20638 degrees of freedom
## Multiple R-squared: 0.4734, Adjusted R-squared: 0.4734
## F-statistic: 1.856e+04 on 1 and 20638 DF, p-value: < 2.2e-16# Nube de puntos con la recta ajustada encima. Pongo los puntos semitransparentes porque
# son más de 20.000 y, si no, se tapan unos a otros.
ggplot(datos, aes(MedInc, MedHouseVal)) +
geom_point(alpha = 0.15, color = "steelblue", size = 0.5) +
geom_smooth(method = "lm", se = FALSE, color = "firebrick") +
scale_y_continuous(labels = dollar) +
labs(title = "Valor de vivienda vs. ingreso medio del distrito",
x = "Ingreso medio (decenas de miles USD)",
y = "Valor de vivienda (USD)") +
theme_minimal()
# Guardo intercepto y pendiente aparte para usarlos en la ecuación y la interpretación.
b0 <- coef(modelo_simple)[1]
b1 <- coef(modelo_simple)[2]La recta estimada queda así: valor de la vivienda ≈ 4.5086^{4} + 4.1794^{4} × ingreso medio.
El intercepto ($45,085.58) sería el valor esperado de una vivienda en un distrito con ingreso cero; no tiene lectura práctica, es solo el punto donde arranca la recta. Lo interesante es la pendiente: por cada unidad más de ingreso medio —cada $10,000 adicionales— el valor de la vivienda sube en promedio $41,793.85.
¿Por qué esos números y no otros? Porque el método de mínimos cuadrados elige el intercepto y la pendiente de modo que la suma de los residuos al cuadrado sea la menor posible; dicho de otra forma, busca la recta que deja la menor distancia vertical total entre lo observado y lo que predice. Esa cantidad que se minimiza —la suma de cuadrados de los residuos— es simplemente la suma de las diferencias entre cada valor observado y su predicho, elevadas al cuadrado.
Para saber si el ingreso realmente influye, contrastamos la hipótesis de que la pendiente vale cero (el ingreso no aportaría nada) frente a la de que es distinta de cero.
# tidy() ordena la salida del modelo en una tabla: estimación, error estándar, estadístico
# t y valor-p de cada coeficiente. El valor-p del ingreso es la prueba de significancia.
tidy(modelo_simple)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 45085.58 | 1322.8716 | 34.0816 | 0 |
| MedInc | 41793.85 | 306.8058 | 136.2225 | 0 |
El valor-p del ingreso es prácticamente cero (< 2e-16), muy por debajo del 5%, así que rechazamos la hipótesis nula: el ingreso medio tiene un efecto estadísticamente significativo sobre el valor de la vivienda. La conclusión es clara y va en la dirección esperada —los distritos con mayores ingresos tienen viviendas más caras—.
# R2: proporción de la variabilidad del precio que el modelo logra explicar.
r2_simple <- summary(modelo_simple)$r.squared
r2_simple## [1] 0.4734475El R² vale 0.473, lo que significa que el ingreso por sí solo explica cerca del 47% de la variación en el precio de las viviendas. Es un poder explicativo moderado: razonable para un único predictor —responde de forma decente—, pero deja claro que falta información y que un modelo con más variables debería mejorarlo.
5 Verificación de supuestos (modelo simple)
Que el modelo ajuste no basta: para que la inferencia (valores-p, intervalos) sea válida, el modelo lineal visto en clase exige cuatro supuestos sobre el error: que la relación sea lineal, que los errores sean independientes, que tengan varianza constante (homocedasticidad) y que se distribuyan de forma normal. Los reviso uno por uno apoyándome en los cuatro gráficos de diagnóstico y en algunas pruebas formales.
# Los cuatro gráficos estándar de un lm: residuos vs ajustados (linealidad), Q-Q
# (normalidad), escala-localización (varianza constante) y residuos vs leverage (puntos
# influyentes). Los pongo en una cuadrícula 2x2.
par(mfrow = c(2, 2))
plot(modelo_simple)
par(mfrow = c(1, 1))Linealidad. En el gráfico de residuos contra valores ajustados, la línea de tendencia debería ser plana alrededor de cero. Aquí la nube está bien centrada, pero aparece una banda diagonal arriba a la derecha: son las 965 viviendas topadas en $500,001. Salvo por ese tope, la relación es razonablemente lineal, lo que respalda usar una recta como primera aproximación.
Independencia. El supuesto pide que los errores no estén correlacionados entre sí; lo evalúo con la prueba de Durbin-Watson.
# Durbin-Watson: estadístico entre 0 y 4; cerca de 2 indica errores independientes.
dwtest(modelo_simple)##
## Durbin-Watson test
##
## data: modelo_simple
## DW = 0.65453, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Un estadístico cercano a 2 indicaría independencia. Vale la pena recordar que estos son datos espaciales —distritos vecinos se parecen—, así que es esperable algo de autocorrelación; el supuesto se cumple solo de manera aproximada y lo dejo anotado como limitación, no como un error del modelo.
Homocedasticidad. Contrasto con Breusch-Pagan si la varianza del error es constante (esa es la hipótesis nula).
# Breusch-Pagan: H0 = varianza constante. Un valor-p pequeño la rechaza (heterocedasticidad).
bptest(modelo_simple)##
## studentized Breusch-Pagan test
##
## data: modelo_simple
## BP = 138.86, df = 1, p-value < 2.2e-16El valor-p es prácticamente cero, así que rechazamos la varianza constante: hay heterocedasticidad, algo que también se ve en el gráfico de escala-localización (la dispersión de los residuos crece con el valor ajustado). Esto no sesga los coeficientes, pero sí invalida los errores estándar y los valores-p clásicos, por lo que una inferencia fina debería apoyarse en errores estándar robustos.
Normalidad de los residuos. Aquí hay que elegir bien la herramienta. Con 20,640 observaciones, las pruebas de hipótesis de normalidad (Shapiro —limitada a 5000 datos—, Lilliefors, Anderson-Darling, Jarque-Bera) rechazan la normalidad casi siempre, porque su potencia crece con el tamaño de muestra y terminan detectando desviaciones irrelevantes: el p-valor deja de ser informativo. El diagnóstico robusto no se basa en ese p-valor sino en los coeficientes de forma —la asimetría y la curtosis—, que miden la magnitud de la desviación y no su significancia. Para una normal la asimetría vale 0 y la curtosis 3.
# Coeficientes de forma de los residuos. Son la vía robusta para juzgar normalidad cuando
# n es grande: describen la magnitud de la desviación, no si es "significativa".
res <- residuals(modelo_simple)
c(asimetria = skewness(res),
curtosis = kurtosis(res),
exceso_curtosis = kurtosis(res) - 3)## asimetria curtosis exceso_curtosis
## 1.191219 5.259881 2.259881# Histograma de los residuos con la curva normal teórica superpuesta, como apoyo visual
# para confirmar lo que dicen los coeficientes de forma.
ggplot(tibble(res), aes(res)) +
geom_histogram(aes(y = after_stat(density)), bins = 60,
fill = "steelblue", alpha = 0.6) +
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)),
color = "firebrick", linewidth = 1) +
labs(title = "Residuos vs. curva normal teórica",
x = "Residuo (USD)", y = "Densidad") +
theme_minimal()
Como regla práctica, una asimetría menor a 0.5 en valor absoluto indica simetría y un exceso de curtosis cercano a cero indica colas normales. Lo que vemos se aleja de eso: la asimetría ronda 1.19 (cola larga hacia la derecha, otra vez por las viviendas caras y el tope de $500,001) y la curtosis ronda 5.3 (exceso de unos 2.3), es decir, una distribución leptocúrtica con colas más pesadas que la normal. El Q-Q plot lo confirma: ajusta bien en el centro y se desvía en los extremos. Los residuos, entonces, no son estrictamente normales, pero la desviación es de una forma conocida y moderada. Gracias al Teorema del Límite Central, con esta muestra tan grande la inferencia sobre los coeficientes sigue siendo válida; lo que conviene tomar con cautela son los intervalos de predicción, que sí dependen de que el error sea normal.
En síntesis: la linealidad y la normalidad se cumplen de forma razonable, mientras que la heterocedasticidad es el supuesto que más claramente se rompe.
6 Regresión lineal múltiple
El paso siguiente es un modelo que sume al ingreso las otras tres variables, para ver si explican el precio mejor que el ingreso solo. Como entre los predictores casi no hay correlación, no esperamos problemas de multicolinealidad al meterlos todos juntos.
# Regresión múltiple con los cuatro predictores. Cada coeficiente mide el efecto de su
# variable sobre el precio manteniendo constantes las demás (efecto "parcial").
modelo_multiple <- lm(MedHouseVal ~ MedInc + HouseAge + AveRooms + AveOccup,
data = datos)
summary(modelo_multiple)##
## Call:
## lm(formula = MedHouseVal ~ MedInc + HouseAge + AveRooms + AveOccup,
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -600529 -53483 -14912 36482 747373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3138.89 2198.71 1.428 0.153
## MedInc 44334.63 312.94 141.672 <2e-16 ***
## HouseAge 1687.52 45.17 37.359 <2e-16 ***
## AveRooms -2734.92 241.39 -11.330 <2e-16 ***
## AveOccup -446.06 53.96 -8.267 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 80480 on 20635 degrees of freedom
## Multiple R-squared: 0.5137, Adjusted R-squared: 0.5136
## F-statistic: 5450 on 4 and 20635 DF, p-value: < 2.2e-16# Guardo los coeficientes para escribir la ecuación e interpretarlos en dólares.
cm <- coef(modelo_multiple)La ecuación estimada queda: valor de la vivienda ≈ 3139 + 4.4335^{4} × ingreso + 1688 × antigüedad − 2735 × habitaciones por hogar − 446 × ocupantes por hogar.
Manteniendo constantes las demás variables, cada coeficiente se lee así: cada $10,000 más de ingreso medio suben el valor de la vivienda en $44,334.63; cada año adicional de antigüedad lo aumenta en $1,687.52; cada habitación promedio de más por hogar lo reduce en $2,734.92; y cada persona adicional por hogar lo baja en $446.06. Los dos efectos negativos tienen sentido: más habitaciones por hogar suele asociarse a zonas rurales más baratas, y más ocupantes por vivienda apunta a hacinamiento, típico de distritos con precios bajos.
Estos números no se pueden comparar directamente entre sí, porque cada variable está en su propia escala (dólares, años, conteos por hogar). Para ver cuál pesa más, estandarizo todas las variables y vuelvo a estimar:
# Reajusto el modelo con todo estandarizado (media 0, desviación 1). Así los coeficientes
# quedan en la misma unidad y se pueden comparar para ver cuál tiene mayor efecto.
modelo_std <- lm(scale(MedHouseVal) ~ scale(MedInc) + scale(HouseAge) +
scale(AveRooms) + scale(AveOccup), data = datos)
sort(abs(coef(modelo_std)[-1]), decreasing = TRUE)## scale(MedInc) scale(HouseAge) scale(AveRooms) scale(AveOccup)
## 0.72990540 0.18404856 0.05863877 0.04014738El coeficiente estandarizado más grande vuelve a ser el del ingreso medio, así que es la variable con mayor efecto sobre el precio de las viviendas, por encima de las otras tres.
# Reviso la significancia de cada coeficiente: marco con TRUE los que tienen valor-p < 0.05.
tidy(modelo_multiple) %>%
mutate(significativo = p.value < 0.05)| term | estimate | std.error | statistic | p.value | significativo |
|---|---|---|---|---|---|
| (Intercept) | 3138.8902 | 2198.71406 | 1.427603 | 0.1534214 | FALSE |
| MedInc | 44334.6251 | 312.93743 | 141.672489 | 0.0000000 | TRUE |
| HouseAge | 1687.5214 | 45.17038 | 37.359022 | 0.0000000 | TRUE |
| AveRooms | -2734.9165 | 241.38744 | -11.329987 | 0.0000000 | TRUE |
| AveOccup | -446.0629 | 53.95516 | -8.267289 | 0.0000000 | TRUE |
Los cuatro coeficientes tienen valor-p por debajo de 0.05: todos son significativos al 5%, ninguno sobra.
# Comparo el R2 del modelo simple con el R2 y el R2 ajustado del múltiple. El ajustado
# penaliza añadir variables, así que es el indicado para comparar modelos de distinto tamaño.
r2_mult <- summary(modelo_multiple)$r.squared
r2_adj_mult <- summary(modelo_multiple)$adj.r.squared
tibble(
modelo = c("Simple (MedInc)", "Múltiple"),
R2 = c(r2_simple, r2_mult),
R2_aj = c(NA, r2_adj_mult)
)| modelo | R2 | R2_aj |
|---|---|---|
| Simple (MedInc) | 0.4734475 | NA |
| Múltiple | 0.5137126 | 0.5136183 |
El modelo múltiple alcanza un R² de 0.514 y un R² ajustado de 0.514, frente al 0.473 del modelo simple. Explica algo más de la variabilidad del precio, así que el múltiple es el que mejor responde, aunque la mejora es moderada: el ingreso sigue siendo el factor dominante. Comparo con el R² ajustado a propósito, porque penaliza meter variables de más; como sube respecto al simple, las tres variables nuevas aportan información real y no solo inflan el ajuste.
7 Verificación de supuestos (modelo múltiple)
Repito el diagnóstico para el modelo múltiple y agrego el supuesto que solo aparece cuando hay varios predictores: que no exista multicolinealidad entre ellos.
# Mismos cuatro gráficos de diagnóstico, ahora para el modelo múltiple.
par(mfrow = c(2, 2))
plot(modelo_multiple)
par(mfrow = c(1, 1))El VIF mide cuánto se infla la varianza de cada coeficiente por culpa de su correlación con los demás predictores. La regla práctica vista en clase dice que un VIF por encima de 10 (incluso de 5) ya es señal de problema.
# Factor de inflación de la varianza por predictor. Cerca de 1 = sin multicolinealidad.
vif(modelo_multiple)## MedInc HouseAge AveRooms AveOccup
## 1.126350 1.029877 1.136638 1.000690Todos los VIF están cerca de 1, así que no hay multicolinealidad. Esto confirma lo que ya anticipaban las correlaciones (la mayor entre predictores era apenas 0.327) y nos deja interpretar cada coeficiente por separado con tranquilidad.
# Reviso los otros tres supuestos de una vez: homocedasticidad (Breusch-Pagan),
# independencia (Durbin-Watson) y normalidad (coeficientes de forma de los residuos).
bptest(modelo_multiple)##
## studentized Breusch-Pagan test
##
## data: modelo_multiple
## BP = 644.69, df = 4, p-value < 2.2e-16dwtest(modelo_multiple)##
## Durbin-Watson test
##
## data: modelo_multiple
## DW = 0.80699, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0res_m <- residuals(modelo_multiple)
c(asimetria = skewness(res_m), curtosis = kurtosis(res_m))## asimetria curtosis
## 1.157657 5.862206El panorama del resto de supuestos es el mismo que en el modelo simple: Breusch-Pagan sigue rechazando la varianza constante (la heterocedasticidad persiste), Durbin-Watson refleja la dependencia espacial, y los coeficientes de forma (asimetría ≈ 1.16, curtosis ≈ 5.9) confirman residuos sesgados a la derecha y con colas pesadas. En otras palabras, sumar variables mejora el ajuste pero no corrige las violaciones de supuestos, que vienen de la naturaleza de los datos y no de cómo planteamos el modelo.
Y aquí vale la pena cerrar el círculo con lo que vimos en el EDA: buena parte de esa falta de normalidad y de varianza constante apunta directamente a los datos extremos. Los distritos con ocupación o habitaciones por hogar disparatadas, junto con el tope de $500,001 en el precio, son justo los que alargan las colas de los residuos y empujan su varianza hacia arriba. El modelo no falla por estar mal planteado, sino porque esos pocos puntos con comportamiento extraño arrastran el diagnóstico.
¿Qué hacer ante esto? La salida más efectiva sería tratar primero esos valores extremos, de forma coherente con el contexto, como recomendé antes (recortarlos o filtrar los imposibles) y reajustar el modelo. De forma complementaria, como la heterocedasticidad afecta la inferencia pero no sesga los coeficientes, se pueden reportar errores estándar robustos para que los valores-p sigan siendo válidos, o transformar el precio con un logaritmo para estabilizar la varianza. Para este taller mantengo el modelo en su escala original, porque permite interpretar los coeficientes directamente en dólares, y dejo anotadas estas limitaciones.
8 Predicción
Por último, estimo el valor de una vivienda en un distrito con ingreso medio de 5, viviendas de 20 años de antigüedad, 6 habitaciones por hogar y 3 ocupantes por hogar.
# Predigo para el distrito pedido y pido dos intervalos al 95%:
# - "confidence": rango para el valor MEDIO de distritos con esas características.
# - "prediction": rango para UN distrito concreto; es más ancho porque suma la
# variabilidad individual, no solo el error de la estimación.
nuevo <- tibble(MedInc = 5, HouseAge = 20, AveRooms = 6, AveOccup = 3)
pred_media <- predict(modelo_multiple, nuevo,
interval = "confidence", level = 0.95)
pred_indiv <- predict(modelo_multiple, nuevo,
interval = "prediction", level = 0.95)
pred_media## fit lwr upr
## 1 240814.8 239366.7 242262.8pred_indiv## fit lwr upr
## 1 240814.8 83064.61 398564.9El valor estimado para ese distrito es de $240,815. Con un 95% de confianza, el valor medio de las viviendas de distritos parecidos cae entre $239,367 y $242,263. Si en cambio queremos el rango para un distrito concreto, el intervalo es bastante más ancho ($83,064.61 a $398,565), porque además del error de la estimación incorpora la variabilidad propia de cada distrito.
9 Conclusiones
El ingreso medio es, de lejos, la variable más ligada al valor de las viviendas (r = 0.688) y el predictor con mayor efecto tanto en el modelo simple como en el múltiple. El modelo simple ya captura buena parte de la señal (R² = 0.47) y el múltiple la mejora hasta R² = 0.51 con sus cuatro predictores, todos significativos al 5%; por eso el múltiple es el que mejor explica el precio, aunque la diferencia no sea enorme.
El diagnóstico, eso sí, deja dos advertencias: hay heterocedasticidad (la varianza del error crece con el precio) y los datos están censurados en el tope de $500,001, lo que de paso explica las colas pesadas de los residuos. Para el distrito propuesto —ingreso 5, antigüedad 20, 6 habitaciones y 3 ocupantes por hogar— el valor estimado es de $240,815.
De cara a mejorar el trabajo, los pasos naturales serían tres, en orden de utilidad: transformar el precio con un logaritmo (atacaría a la vez la heterocedasticidad y la asimetría de los residuos), tratar explícitamente los valores topados —por ejemplo excluyéndolos o modelándolos como censura—, y usar errores estándar robustos para que la inferencia no dependa del supuesto de varianza constante. Con eso el modelo ganaría tanto en validez de los supuestos como en capacidad explicativa. ```