Predicción de precios de las casas en Boston

Autor/a

Pedro Albarrán

Introducción

Comentario General

Este documento debe entenderse como un ejemplo, no la guía única para el proyecto.

Aquí muestro algunas cosas que en vuestros trabajos probablemente NO queréis mostrar. Por ejemplo, he incluido la opción de que se muestre o se oculte el código de todo lo que he hecho, pero vosotros debéis pensar qué y cuándo queréis mostrar algo. También muestro algunos resultados que, como se ha discutido en clase, probablemente tampoco queráis mostrar (o no de la misma manera).
He omitido muchos detalles de algunas fases del trabajo que ya se han ido comentando y trabajando en clase (características de los datos, procedimientos, tablas con encabezados adecuados, gráficos con ejes correctamente nombrados, comentario de resultados, etc.)
Cada conjunto de datos es diferente y cada análisis es diferente. Se requiere un tratamiento distinto de los datos: (eliminar o imputar valores ausentes, transformar variables, agrupar categorías, discretizar, etc.), diferentes algoritmos y especificaciones (combinaciones de variables a incluir y sus transformaciones). De hecho, se espera que probéis distintas opciones y discutáis los resultados.

En resumen, este ejemplo incluye pruebas, gráficos y tablas que vosotros no incluiréis. También deberéis probar otras que aquí no he probado. Pero también se espera discusiones algo más detalladas (sin ser excesivas).

Introducción y Objetivos

En este trabajo, se analizará un conjunto de datos con información sobre precios y otros atributos de una muestra de viviendas en Boston. Por un lado, el objetivo es examinar la influencia de varios atributos del vecindario en los precios de la vivienda, en un intento por descubrir las variables explicativas más adecuadas. Por otro lado, la construcción de un modelo de predicción permitirá determinar el valor por el que se puede poner en el mercado una vivienda o detectar si alguna está infravalorada o sobrevalorada dadas sus características.

El objetivo del trabajo es un elemento importante que debéis desarrollar adecuadamente. Todo el análisis y comentarios posteriores deben estar orientados a responder a las preguntas planteados por los objetivos.

Datos

En este trabajo vamos a utilizar un conjunto de datos, “The Boston Housing Price”, derivados de la información recopilada por el Servicio de Censos de los Estados Unidos sobre las viviendas en el área de Boston (Massachusetts). Podemos encontrar algunos detalles adicionales sobre estos datos aquí. Los datos pueden obtenerse desde esa misma página, pero por sencillez los leemos en formato de valores separados por comas (CSV) desde aquí.

Boston <- read_csv("data/BostonHousing.csv")

Los datos tienen 506 observaciones y 14 variables. (NOTA: NO está claro que queráis mostrar este código.)

dim(Boston)

[1] 506 14

Una descripción completa los atributos disponibles puede encontrarse en Apéndice A. Estos datos fueron originalmente utilizados en un estudio sobre el impacto de la contaminación del aire (utilizando las concentraciones de óxido de nitrógeno). En este trabajo, consideramos el efecto de otras características de la zona donde se encuentra la casa como la proximidad al río Charles, la distancia a los principales centros de empleo, la calidad de las escuelas (medida por la ratio de número de alumnos por maestro) y los niveles de delincuencia. Nuestra variable de interés para predecir es el valor mediano del precio de la vivienda en mil dólares (denotado por MEDV). Notad que nos centramos en las viviendas ocupadas por sus propietarios, es decir, consideramos que el valor de las casas destinadas al alquiler se determina según un proceso diferente.

Exploración Inicial

Debemos considerar si los datos están listos para trabajar o requieren algún tipo de limpieza, ordenación o transformación. En primer lugar, comprobamos si el tipo de datos de cada variable es el adecuado. En este caso, todas las variables son numéricas, lo cual se corresponde con la información cuantitativa de la mayoría de ellas. La variable que nos dice si la zona de la casa está cerca del río es una variable binaria, es decir, aporta información cualitativa. En este caso, no es crucial convertirla en un factor (en los modelos y para otras cuestiones, convertimos los factores en variables binarias). En otros casos, podemos necesitar convertir más variables con información cualitativa a factores o eliminar variables de tipo carácter o con información que no podamos procesar. NOTA: para esto habremos comprobado el tipo de cada variable, habremos mirado los valores y contrastado con la información que según su descripción debería tener. Utilizaríamos un código como el siguiente, aunque probablemente NO queráis incluir en el documento ni lo uno ni lo otro.

glimpse(Boston)

Rows: 506
Columns: 14
$ crim    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829,…
$ zn      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 1…
$ indus   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
$ chas    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ nox     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
$ rm      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631,…
$ age     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, 9…
$ dis     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.9505…
$ rad     <dbl> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,…
$ tax     <dbl> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 31…
$ ptratio <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
$ b       <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.90…
$ lstat   <dbl> 4.98, 9.14, 4.03, 2.94, 5.33, 5.21, 12.43, 19.15, 29.93, 17.10…
$ medv    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15…

head(Boston) %>% kbl() %>% kable_paper("hover")

crim	zn	indus	nox	rm	age	dis	rad	tax	ptratio	b	lstat	medv
0.00632	18	2.31	0.538	6.575	65.2	4.0900	1	296	15.3	396.90	4.98	24.0
0.02731	0	7.07	0.469	6.421	78.9	4.9671	2	242	17.8	396.90	9.14	21.6
0.02729	0	7.07	0.469	7.185	61.1	4.9671	2	242	17.8	392.83	4.03	34.7
0.03237	0	2.18	0.458	6.998	45.8	6.0622	3	222	18.7	394.63	2.94	33.4
0.06905	0	2.18	0.458	7.147	54.2	6.0622	3	222	18.7	396.90	5.33	36.2
0.02985	0	2.18	0.458	6.430	58.7	6.0622	3	222	18.7	394.12	5.21	28.7

Boston <- Boston %>% mutate(chas = factor(chas, 
                         levels=c(0,1), labels = c("No", "Yes")))

(NOTA: los nombres de las columnas y de las filas/variables son “mejorables”: si el documento está en castellano las columnas no deberían tener nombres en inglés y sería preferible que apareciera un nombre más descriptivo de las variables. Pero no voy a exigir esto en el plazo que tenemos. Eso sí, al menos que en algún sitio aparezca la descripción completa del nombre abreviado de las variables, como hago aquí con el Apéndice A.)

Análisis Exploratorio y Visualización de los datos

La siguiente fase consiste en analizar la distribución de valores de cada variable (análisis de variación) y las posibles relaciones entre ellas (análisis de covariación). Esto nos puede llevar a realizar limpieza adicional de los datos (en particular, relacionada con valores ausentes y puede que atípicos) o transformaciones de los datos (como tomar logaritmos o discretizar alguna variable). También podemos encontrar características de los datos que sean de interés por sí mismas como para especificar los modelos. Para este análisis nos podemos ayudar en bibliotecas que realizan algunas de las tareas de forma automatizada. PERO recordad que NO queremos en general mostrar la salida directa de estos paquetes, sino que la utilizaremos para aprender nosotros y luego mostrar aquello que consideremos más interesante.

Por ejemplo, usando la biblioteca skimr podemos ver

skim(Boston)

Data summary
Name	Boston
Number of rows	506
Number of columns	14
_______________________
Column type frequency:
factor	1
numeric	13
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
chas	0	1	FALSE	2	No: 471, Yes: 35

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
crim	1	3.61	8.60	0.01	0.08	0.26	3.68	88.98	▇▁▁▁▁
zn	1	11.36	23.32	0.00	0.00	0.00	12.50	100.00	▇▁▁▁▁
indus	1	11.14	6.86	0.46	5.19	9.69	18.10	27.74	▇▆▁▇▁
nox	1	0.55	0.12	0.38	0.45	0.54	0.62	0.87	▇▇▆▅▁
rm	1	6.28	0.70	3.56	5.89	6.21	6.62	8.78	▁▂▇▂▁
age	1	68.57	28.15	2.90	45.02	77.50	94.07	100.00	▂▂▂▃▇
dis	1	3.80	2.11	1.13	2.10	3.21	5.19	12.13	▇▅▂▁▁
rad	1	9.55	8.71	1.00	4.00	5.00	24.00	24.00	▇▂▁▁▃
tax	1	408.24	168.54	187.00	279.00	330.00	666.00	711.00	▇▇▃▁▇
ptratio	1	18.46	2.16	12.60	17.40	19.05	20.20	22.00	▁▃▅▅▇
b	1	356.67	91.29	0.32	375.38	391.44	396.22	396.90	▁▁▁▁▇
lstat	1	12.65	7.14	1.73	6.95	11.36	16.96	37.97	▇▇▅▂▁
medv	1	22.53	9.20	5.00	17.02	21.20	25.00	50.00	▂▇▅▁▁

o usando dlookr

Boston %>% 
  describe() %>%
  select(described_variables:kurtosis) %>% 
  kbl(digits = 2) %>% kable_paper("hover")

described_variables	n	mean	sd	se_mean	IQR	skewness	kurtosis
crim	506	3.61	8.60	0.38	3.60	5.22	37.13
zn	506	11.36	23.32	1.04	12.50	2.23	4.03
indus	506	11.14	6.86	0.30	12.91	0.30	-1.23
nox	506	0.55	0.12	0.01	0.17	0.73	-0.06
rm	506	6.28	0.70	0.03	0.74	0.40	1.89
age	506	68.57	28.15	1.25	49.05	-0.60	-0.97
dis	506	3.80	2.11	0.09	3.09	1.01	0.49
rad	506	9.55	8.71	0.39	20.00	1.00	-0.87
tax	506	408.24	168.54	7.49	387.00	0.67	-1.14
ptratio	506	18.46	2.16	0.10	2.80	-0.80	-0.29
b	506	356.67	91.29	4.06	20.85	-2.89	7.23
lstat	506	12.65	7.14	0.32	10.01	0.91	0.49
medv	506	22.53	9.20	0.41	7.98	1.11	1.50

Boston %>% 
  describe() %>%
  select(described_variables, p00:p100) %>% 
  kbl(digits = 2) %>% kable_paper("hover")

described_variables	p00	p01	p05	p10	p20	p25	p30	p40	p50	p60	p70	p75	p80	p90	p95	p99	p100
crim	0.01	0.01	0.03	0.04	0.06	0.08	0.10	0.15	0.26	0.55	1.73	3.68	5.58	10.75	15.79	41.37	88.98
zn	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	12.50	20.00	42.50	80.00	90.00	100.00
indus	0.46	1.25	2.18	2.91	4.39	5.19	5.96	7.38	9.69	12.83	18.10	18.10	18.10	19.58	21.89	25.65	27.74
nox	0.38	0.40	0.41	0.43	0.44	0.45	0.47	0.51	0.54	0.58	0.60	0.62	0.67	0.71	0.74	0.87	0.87
rm	3.56	4.52	5.31	5.59	5.84	5.89	5.95	6.09	6.21	6.38	6.50	6.62	6.75	7.15	7.59	8.33	8.78
age	2.90	6.61	17.72	26.95	37.80	45.02	52.40	65.40	77.50	85.90	91.80	94.07	95.60	98.80	100.00	100.00	100.00
dis	1.13	1.21	1.46	1.63	1.95	2.10	2.26	2.64	3.21	3.88	4.54	5.19	5.62	6.82	7.83	9.22	12.13
rad	1.00	1.00	2.00	3.00	4.00	4.00	4.00	5.00	5.00	5.00	8.00	24.00	24.00	24.00	24.00	24.00	24.00
tax	187.00	188.00	222.00	233.00	273.00	279.00	289.00	307.00	330.00	398.00	437.00	666.00	666.00	666.00	666.00	666.00	711.00
ptratio	12.60	13.00	14.70	14.75	16.60	17.40	17.80	18.40	19.05	19.70	20.20	20.20	20.20	20.90	21.00	21.20	22.00
b	0.32	6.73	84.59	290.27	364.31	375.38	378.66	387.97	391.44	393.53	395.47	396.22	396.90	396.90	396.90	396.90	396.90
lstat	1.73	2.88	3.71	4.68	6.29	6.95	7.77	9.53	11.36	13.33	15.62	16.96	18.06	23.04	26.81	33.92	37.97
medv	5.00	7.01	10.20	12.75	15.30	17.02	18.20	19.70	21.20	22.70	24.15	25.00	28.20	34.80	43.40	50.00	50.00

En este caso, el conjunto de datos está ya bastante limpio, por lo que apenas necesitamos realizar cambios y podemos dedicar más tiempo al resto del proceso. En otros casos, deberemos realizar más trabajo en esta parte; en particular, relacionado con valores ausentes, errores en los datos, etc. como hemos visto a lo largo del curso y, en particular, en el tema de análisis exploratorio de datos.

Análisis de variación

Como primer elemento a destacar, estos datos no contiene valores ausentes en ninguna de las variables. En caso contrario, deberíamos identificar cuántas observaciones y qué variables están afectadas. Sabemos que podemos posponer la imputación de valores a una fase posterior (como un paso del pre-procesado antes de estimar un modelo), pero es conveniente tener una visión general y pensar si algunas observaciones probablemente serán descartadas (si tienen muchos valores ausentes y sobre todo afectan a la variable dependiente).

Podemos centrarnos en describir con más detalles algunas distribuciones. Esto nuevamente es un EJEMPLO dependiendo de las variables que tengamos y de qué observemos. En general, caracterizar la variable dependiente suele ser una buena idea. Visualizamos la distribución y densidad del precio mediano de las viviendas. La curva negra representa la densidad. Vemos que el valor medio del precio de la vivienda está sesgado a la derecha. Es decir, observamos precios muy altos con una frecuencia mayor de la esperada en una distribución simétrica donde existiría la misma proporción por encima y debajo de la media.

Boston %>%  ggplot(aes(x=medv)) + geom_histogram(aes(y=..density..))+ geom_density() + ggtitle("Distribución del Precio") + xlab("Precio de las casas") + ylab("Densidad")

Figura 1. Distribución del precio de la vivienda

Dada esta asimetría, quizás debamos considerar modelizar posteriormente esta variable transformada en logaritmos. La razón: se aprecia un comportamiento que puede modelizarse mejor de forma no lineal. También se puede nota una acumulación de valores en 50 mil dólares. Se puede observar en los resultados de describe() que ese valor exacto se repite varias veces, NO es producto de la discretización del gráfico en la que se acumulan varios valores diferentes en un intervalo en torno a 50; también debemos probar distintos anchos de intervalo como se ha discutido en clase.

También podemos representar gráficamente o en un tabla la única variable categórica que tenemos. La conclusión no es particularmente interesante: solo unas pocas zonas de la ciudad están cerca del río.

Boston %>%  count(chas) %>% 
  mutate(freq=n/sum(n)) %>% 
  kbl(col.names=c("Casa cercana al río","Número de casos", "Frecuencia"),
                  caption = "Tabla 1. Distribución de Casas según cercanía al río") %>% kable_paper("hover")

Tabla 1. Distribución de Casas según cercanía al río
Casa cercana al río	Número de casos	Frecuencia
No	471	0.93083
Yes	35	0.06917

En el caso de variables binarias las podemos representar de varias maneras: como una distribución o con una sola barra (NOTA: los gráficos siguientes son redundantes en esta caso, con uno de ellos sería más que suficiente en caso de considerar relevante esta información.)

Boston %>%  ggplot() + geom_bar(aes(x=chas)) +  xlab("Zona cercana al río") +ylab("Número de casos")

Boston %>%  ggplot() + geom_bar(aes(x="",fill=chas)) + labs(fill="Zona cercana al río") +ylab("Número de casos")

Figura 2. Distribución de la cercanía al río

También podemos mostrar algunas otras características interesantes mediante gráficos y/o tablas de estadísticos descriptivos. Algunas variables como el número de habitaciones tienen distribuciones bastante simétricas. Mientras que otras, como la edad o el porcentaje de población desfavorecida muestran claras asimetrías: hay una alta concentración de casas “viejas” y de zonas no desfavorecidas. NOTA: Recordad que habría que probar varios anchos de intervalos (binwidth) para asegurarnos de entender la forma de la distribución. También debéis poner nombres suficientemente descriptivos e informativos a los gráficos, los ejes, la leyenda, etc. (Quizás no es el caso en algunos de los que presento aquí).

Boston %>%  ggplot(aes(x=age)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Edad") + ylab("Densidad")

Boston %>%  ggplot(aes(x=lstat)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Porcentaje de población desfavorecida") + ylab("Densidad")

Notad que en la variable edad nuevamente hay una concentración de valores en 100; en la salida describe() mostrada anteriormente, se aprecia mejor que ese valor exacto está en los datos originales, no resulta de que se agrupen valores en el gráfico.

El caso de la distancia a los centros de empleo es similar a las dos anteriores: una gran concentración en zonas bien conectadas, aunque una cola de zonas alejadas. Se podría omitir: no hay que mostrar gráficos o tablas de cada variable ni comentar necesariamente las características de la distribución de todas, solo de aquellas con rasgos interesante o relevantes.

Otra variable en principio relacionada, el índice de accesibilidad, muestra una distribución “poco continua”: además de un cúmulo de valores en la cola derecha, hay muchos huecos vacíos. Si probáis un transformación logarítmica, veréis que no cambia en esencia. Las variables con este forma en su distribución suelen ser candidatas a ser discretizadas.

Boston %>%  ggplot(aes(x=dis)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Distancia a centro de trabajo") + ylab("Densidad")

Boston %>%  ggplot(aes(x=rad)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Índice de accesibilidad") + ylab("Densidad")

Algunas variables tienen distribuciones con características poco reseñables: unas con valores distribuidos de forma relativamente homogénea, otras dispersas, con concentraciones en valores aislados en medio o en los extremos de la distribución, pero no aportan mucho información (se podrían omitir). En este caso, quizás se podría notar una concentración de zonas con altos impuestos, muy diferenciadas del resto.

Boston %>%  ggplot(aes(x=nox)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Concentración de óxidos nítricos") + ylab("Densidad")

Boston %>%  ggplot(aes(x=ptratio)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Ratio de alumnos por profesor") + ylab("Densidad")

Boston %>%  ggplot(aes(x=tax)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Impuesto de la propiedad") + ylab("Densidad")

Algo más interesantes son algunas variables que muestran polaridad en sus valores o una excesiva acumulación en algunos. Por ejemplo, la criminalidad y el porcentaje de población de color tienen distribuciones muy asimétricas y, en el segundo caso, persiste incluso tras transformar en logaritmos.

Boston %>%  ggplot(aes(x=crim)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Criminalidad") + ylab("Densidad")
Boston %>%  ggplot(aes(x=crim)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Criminalidad") + ylab("Densidad") + scale_x_log10()

Boston %>%  ggplot(aes(x=b)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Población de color") + ylab("Densidad") 
Boston %>%  ggplot(aes(x=b)) + geom_histogram(aes(y=..density..))+ geom_density() + xlab("Población de color") + ylab("Densidad") + scale_x_log10()

Estas variables y algunas otras anteriores son candidatas a ser discretizadas. La criminalidad, por ejemplo, no solo muestra una concentración en unos pocos valores, sino que una vez transformada en logaritmos se aprecian dos grupos diferenciados, como también pasaba con los impuestos. En el caso de la ratio de profesor/alumno también unos valores con gran concentración de frecuencia y muy pocos por encima de este por lo que podrían agruparse juntos. En el caso de la población de color, vemos que a partir del percentil 75, los valores son prácticamente iguales y antes del percentil 10 son mucho menores que en el resto de la distribución. En estos casos de variables con valores concentrados o infrecuentes y con saltos o huecos, no podemos decir que la variable que observamos en nuestra muestra tenga una apariencia de variable continua, aunque en principio lo sea. Por tanto, agrupar y discretizar es una buena opción. En particular, es más fácil identificar el efecto medio sobre el precio de la vivienda de un rango de valores (ej., zonas de baja criminalidad frente a alta) que el efecto de incrementar en un punto la variable (cuando en los datos no observamos valores con ese punto más). En otras palabras, vamos a modelizar efectos flexibles no lineales.

MUY IMPORTANTE:

NO mostréis todos los gráficos o tablas que se os ocurran o solo porque los veáis aquí. ELEGID adecuadamente aquellos que consideréis más relevantes o interesantes y aporten información útil para responder a los objetivos planteados.

Análisis de covariación

Empezamos analizando la relación entre nuestra variable de interés y la única variable categórica que tenemos inicialmente, para lo que podríamos presentar alguna (NO todas) de las siguientes figuras

Boston %>% ggplot(aes(y = medv, x = chas)) +  geom_boxplot() + ylab("Densidad") + xlab("Cerca del río") + ylab("Precio") 
Boston %>% ggplot(aes(x = medv)) + geom_density(mapping = aes(colour = chas)) + xlab("Precio") + ylab("Densidad")  + labs(color = "Cerca del río")

Figura 7a. Distribución del precio por cercanía al rio

Boston %>% ggplot(aes(x = medv)) + geom_density() + facet_wrap(~chas) + ylab("Densidad")  + xlab("Precio")

Figura 7b. Distribución del precio, según cercanía al rio

Parece que las casa cercanas al río tienen un precio superior, aunque la diferencia no parece grande. Ambas distribuciones son asimétricas, aunque en el caso de casas cercanas al río la cola derecha no es tan larga. Mediante una regresión simple o calculando las medias podemos comprobar si existen diferencias en media y si son significativas:

lm(data = Boston, medv ~ chas) %>% broom::tidy() %>%  kbl(digits = 2, caption = "Tabla 2a. Precio según cercanía al río") %>% kable_paper("hover")

Tabla 2a. Precio según cercanía al río
term	estimate	std.error	statistic	p.value
(Intercept)	22.09	0.42	52.9	0
chasYes	6.35	1.59	4.0	0

Boston %>% group_by(chas) %>% describe(medv) %>% select(described_variables:n, mean, se_mean) %>%  kbl(digits = 2, caption = "Tabla 2b. Precio según cercanía al río") %>% kable_paper("hover")

Tabla 2b. Precio según cercanía al río
described_variables	chas	n	mean	se_mean
medv	No	471	22.09	0.41
medv	Yes	35	28.44	2.00

También podríamos analizar si la variable binaria de cercanía al río está relacionada. (NOTA: esto nuevamente es un código de ejemplo para hacer de forma fácil este proceso. NO es necesario que vosotros lo hagáis así.)

vars <- c("crim", "zn", "indus", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "b", "lstat")

tabla <-list()
for (var in vars) {
  formula <- paste0(var," ~ chas")
  tabla[[var]] <- lm(data = Boston, formula) %>% tidy() %>% filter(term=="chasYes") %>% select(estimate, p.value)
}

tabla %>% bind_rows(.id="Variable")  %>%  kbl(digits = 2, caption = "Tabla 2c. Diferencias en variable por cercanía al río") %>% kable_paper("hover")

Tabla 2c. Diferencias en variable por cercanía al río
Variable	estimate	p.value
crim	-1.89	0.21
zn	-3.92	0.34
indus	1.70	0.16
nox	0.04	0.04
rm	0.25	0.04
age	9.59	0.05
dis	-0.82	0.03
rad	-0.25	0.87
tax	-23.61	0.42
ptratio	-1.04	0.01
b	17.54	0.27
lstat	-1.52	0.23

Vemos que las casas cercanas al rio son más antiguas, con más habitaciones, más cercanas al centro de trabajo, con más contaminación y mejores condiciones escolares.

A continuación podemos analizar rápidamente si las variables continuas están relacionadas con nuestra variable de interés, precio de la vivienda, y entre ellas. Lo podemos hacer mediante distintos análisis de correlación, en una tabla (excesivamente larga) o visualmente.

Boston %>% correlate() %>% 
  filter(as.integer(var1) > as.integer(var2)) %>% 
  kbl(digits = 2, caption = "Tabla 3. Correlaciones") %>% kable_paper("hover")

Boston %>% correlate() %>% plot() 

Boston %>% mutate(logmedv=log(medv)) %>% select(-medv) %>% correlate() %>% plot()

Figura 8. Correlaciones entre variables continuas

Hemos considerado la correlación tanto con el precio como con su logaritmo, dado lo discutido anteriormente. Sin embargo, apenas se aprecian diferencias.

Vemos que existe una fuerte correlación (positiva o negativa) entre el precio y varias variables que intuitivamente consideraríamos como importantes. El número de habitaciones tiene la correlación positiva más fuerte con el valor medio del precio de la vivienda, mientras que el porcentaje de la población desfavorecida y el número de alumnos por docente tienen una correlación negativa fuerte. También es evidente que las zonas más industriales y la contaminación están fuertemente correlacionados positivamente entre sí, puesto que los niveles de óxido nítrico tienden a aumentar con el aumento de las industrias. También vemos que las zonas con más población desfavorecida son las más industriales y contaminadas, con casas más antiguas y de menos habitaciones y con escuelas con un mayor ratio de alumnos por profesor. Debemos recordar esto de cara a la especificación de los modelos de regresión lineal.

Sin embargo, esto no considera posibles relaciones no lineales. Para ello vamos a representar varios gráficos de dispersión y un ajuste no lineal. Nuevamente, en vuestro trabajo no mostraréis necesariamente todos estos gráficos sino una selección después de haberlos vistos.

vars <- c("crim", "zn", "indus", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "b", "lstat")

for (v in vars){
  graf <- Boston %>% 
            ggplot(aes_string(x = v, y = "medv")) +
            geom_point() +  geom_smooth() +
            labs(x = v, y = "Precio de las casas ($1000s)")
  print(graf)
}

vars <- c("crim", "zn", "indus", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "b", "lstat")

for (v in vars) {
  migraf <-  Boston %>% 
              ggplot(aes_string(x = v, y = "medv")) +
              geom_point() +  geom_smooth() +
              labs(x = v, y = "Precio de las casas ($1000s)")  +
              scale_y_log10() 
  print(migraf)
}

Figura 9b. Gráficos de dispersión (en escala logaritmica)

En primer lugar, no se aprecian grandes diferencias entre el modelo con el precio sin transformar o en logaritmos. En segundo lugar, sí se aprecia cierta no linealidad en la relación con las variables de edad, número de habitaciones y porcentaje de población desfavorecida. En el resto de relaciones, no están tan claras por la acumulación de valores.

También se podrían haber probado si este análisis de covariación es distinto según distintos valores de una varible categórica. Por ejemplo, visualizar la relación entre precio y tasa de pobreza cuando la zona está cerca del rio y cuando no está cerca.

Podemos probar discretizando algunas de las variables comentadas anteriormente. Por ejemplo, hacemos dos grupos de criminalidad; los umbrales para discretizar no tienen una justificación muy formal: se basan en lo que aproximadamente hemos visto.

  Boston %>% 
      mutate(crim.alta = cut(crim, breaks = c(0,1,Inf), 
                             include.lowest = TRUE,
                             labels = c("Baja","Alta") ) ) %>% 
      ggplot(aes(y = medv, x = crim.alta)) +  geom_boxplot() + ylab("Densidad") + xlab("Criminalidad alta") + ylab("Precio")

Figura 10. Gráficos Criminalidad Discreta

  Boston %>% 
      mutate(crim.alta = cut(crim, breaks = c(0,1,Inf), 
                             include.lowest = TRUE,
                             labels = c("Baja","Alta") ) ) %>% 
  lm(data = ., medv ~ crim.alta) %>% tidy() %>%  kbl(digits = 2, caption = "Tabla 4. Precio según criminalidad") %>% kable_paper("hover")

Tabla 4. Precio según criminalidad
term	estimate	std.error	statistic	p.value
(Intercept)	25.11	0.47	53.91	0
crim.altaAlta	-7.50	0.79	-9.44	0

Tanto el gráfico como la regresión apuntan a un efecto significativo de la criminalidad sobre los precios. En principio deberíamos probar con otros puntos de corte para discretizar, pero por simplicidad utilizaremos este obtenido a partir del análisis exploratorio.

Podemos proceder de manera similar con otras variables. Nuevamente, tanto el número de grupos como los valores de corte no se derivan de forma super rigurosa, sino en base al análisis exploratorio. Debería probarse con otras variantes.

Boston %>% 
  mutate(dis.alta = cut(dis, breaks = c(0,3,Inf), 
                        include.lowest = TRUE,
                        labels = c("Baja","Alta") ) ) %>% 
  lm(data = ., medv ~ dis.alta) %>% tidy() %>%
  kbl(digits = 2, caption = "Tabla 5. Precio según distancia") %>% kable_paper("hover")

Tabla 5. Precio según distancia
term	estimate	std.error	statistic	p.value
(Intercept)	19.73	0.57	34.69	0
dis.altaAlta	5.32	0.78	6.79	0

Boston %>% 
  mutate(rad.alta = cut(dis, breaks = c(0,10,Inf), 
                        include.lowest = TRUE,
                        labels = c("Baja","Alta") ) ) %>% 
  lm(data = ., medv ~ rad.alta) %>% tidy() %>%  
  kbl(digits = 2, caption = "Tabla 6. Precio según accesibilidad") %>% kable_paper("hover")

Tabla 6. Precio según accesibilidad
term	estimate	std.error	statistic	p.value
(Intercept)	22.53	0.41	54.79	0.00
rad.altaAlta	-0.21	4.14	-0.05	0.96

Boston %>% 
  mutate(tax.alta = cut(tax, breaks = c(0,350,500,Inf), 
                        include.lowest = TRUE, 
                        labels = c("Baja","Media","Alta")) ) %>%
  lm(data = ., medv ~ tax.alta) %>% tidy() %>%  
  kbl(digits = 2, caption = "Tabla 7. Precio según impuestos") %>% kable_paper("hover")

Tabla 7. Precio según impuestos
term	estimate	std.error	statistic
(Intercept)	25.88	0.51	51.16
tax.altaMedia	-3.68	0.96	-3.83
tax.altaAlta	-9.60	0.87	-11.06

Boston %>% 
  mutate(black.cat = cut(b, breaks = c(0,100, 395,Inf), 
                         include.lowest = TRUE,
                         labels = c("Baja","Media","Alta")) ) %>%
  lm(data = ., medv ~ black.cat) %>% tidy() %>%  
  kbl(digits = 2, caption = "Tabla 8. Precio según población de color") %>% kable_paper("hover")

Tabla 8. Precio según población de color
term	estimate	std.error	statistic
(Intercept)	12.33	1.58	7.82
black.catMedia	11.63	1.65	7.03
black.catAlta	9.48	1.72	5.52

Boston %>% 
  mutate(black.alta = cut(b, breaks = c(0, 100,Inf), 
                          include.lowest = TRUE,
                          labels = c("Baja","Alta")) ) %>% 
  lm(data = ., medv ~ black.alta) %>% tidy() %>%  
  kbl(digits = 2, caption = "Tabla 8b. Precio según población de color") %>% kable_paper("hover")

Tabla 8b. Precio según población de color
term	estimate	std.error	statistic	p.value
(Intercept)	12.33	1.59	7.77	0
black.altaAlta	10.87	1.64	6.64	0

(NOTA: estas tablas se podrían haber hecho con un bucle. También notad que NO he mantenido la variable discretizada en los datos, aunque podría haberlo hecho.)

IMPORTANTE

La fase de análisis exploratorio puede tener resultados interesantes por sí mismos. Es importante comentarlos y discutirlos adecuadamente, en relación con el objetivo del trabajo.

Pero sobre todo esta fase tiene como objetivo aprender de los datos de cara a la siguiente: modelización. Por tanto, es mucho más importante comentar, aquí o en la fase de modelización cuestiones cómo

qué variables incluir en los modelos: aunque esto no es crucial para algunos algoritmos que seleccionan.

b.cómo las vamos a incorporar: en función de los resultados previos podemos queremos incluir transformaciones no lineales, discretizaciones, agrupando categorías, etc.

Modelos

NOTA IMPORTANTE

En general, hemos usuado habitualmente step_log() para transformar en logaritmos la variable dependiente de los modelos. Si al final queremos predecir la variable dependiente es preferible generar la variable en logaritmos y usarla directamente en los modelos.

Boston <- Boston %>% mutate(lmedv = log(medv))

Muestras de entrenamiento y prueba

En el análsisis exploratorio, hemos visto que varias variables explicativas podrían ser transformadas tomando logaritmos o polinomios, discretizando, etc. Todo esas transformaciones se pueden hacer en el pre-procesado de tidymodels.

A continuación, hacemos la partición de los datos reservando una proporción del 80% como conjuntos de datos de entrenamiento y el 20% restante como prueba.

set.seed(1)
bostonPart <- Boston %>% initial_split(prop = .8)

Modelos de regresión lineal

Se pueden probar muchisimo modelos de regresión lineal. Demasiados, como deberiáis saber.¹ Si nuestro objetivo es predecir, es preferible estimar, como hemos discutivos, un modelo (o modelos) más “complejos” mediante regresión regularizada (LASSO y/o “ridge regression”); por ejemplo, un modelos con un polinomio de orden grande o muchas interacciones en lugar de ir probando modelos cuadráticos o cúbicos o solo unas pocas interacciones. Por contra, si nuestro objetivo es interpretar qué variables afectan y cuánto a la variables dependiente, podemos estimar el modelo de regresión lineal después de haber estimado mediante LASSO para aprender qué variables parecen ser relevantes.

En general, no haremos esto, pero aquí probaremos un modelo de regresión lineal para el precio como variable dependiente y todas las variables restantes como variables independientes; es decir, sin ninguna relación no lineal. Estimamos este modelo aquí como una extensión del análisis exploratorio, pero considerando correlaciones parciales (ceteris paribus): el efecto de cada variable, manteniendo el resto constante. Entrenamos el modelo con el conjunto de datos de entrenamiento. A continuación se muestran todos los coeficientes.

lm1_receta <- training(bostonPart) %>%            
  recipe(medv ~ chas + crim + zn + indus + nox + rm + 
           age + dis + rad + tax + ptratio + b + lstat) 

lm1_modelo <- linear_reg(mode= "regression", penalty = 0) %>%
                    set_engine("lm")

lm1_flujo <- workflow() %>%
  add_recipe(lm1_receta) %>%
  add_model(lm1_modelo)

lm1_flujo_est <- lm1_flujo %>% fit(data = training(bostonPart)) 

lm1_flujo_est %>% extract_fit_parsnip() %>% tidy() %>% 
   kbl(digits = 2, caption = "Tabla 9. Modelo de Regresión") %>% kable_paper("hover")

Tabla 9. Modelo de Regresión
term	estimate	std.error	statistic	p.value
(Intercept)	32.41	6.06	5.35	0.00
chasYes	3.15	0.97	3.24	0.00
crim	-0.09	0.04	-2.25	0.02
zn	0.04	0.02	2.58	0.01
indus	0.03	0.07	0.43	0.67
nox	-14.95	4.42	-3.38	0.00
rm	4.07	0.49	8.35	0.00
age	0.00	0.02	-0.33	0.74
dis	-1.44	0.23	-6.14	0.00
rad	0.32	0.07	4.37	0.00
tax	-0.01	0.00	-3.18	0.00
ptratio	-0.88	0.15	-5.72	0.00
b	0.01	0.00	3.45	0.00
lstat	-0.55	0.06	-9.64	0.00

Los resultados del modelo con todas las variables son similares a los anteriores con solo una variable cada vez, aunque la edad NO es significativa (ceteris paribus).

lm2_receta <- training(bostonPart) %>%            
                recipe(medv ~ chas + crim + zn + indus + nox + rm + 
                    age + dis + rad + tax + ptratio + b + lstat) %>% 
                step_log(medv)

lm2_flujo <- workflow() %>%
  add_recipe(lm2_receta) %>%
  add_model(lm1_modelo)

lm2_flujo_est <- lm2_flujo %>% fit(data = training(bostonPart)) 

lm2_flujo_est %>% extract_fit_parsnip() %>% tidy() %>% 
   kbl(digits = 3, caption = "Tabla 10. Modelo de Regresión (en logaritmos)") %>% kable_paper("hover")

Tabla 10. Modelo de Regresión (en logaritmos)
term	estimate	std.error	statistic	p.value
(Intercept)	4.133	0.237	17.438	0.000
chasYes	0.108	0.038	2.854	0.005
crim	-0.012	0.002	-7.123	0.000
zn	0.001	0.001	1.821	0.069
indus	0.002	0.003	0.785	0.433
nox	-0.714	0.173	-4.135	0.000
rm	0.083	0.019	4.353	0.000
age	0.000	0.001	-0.154	0.877
dis	-0.051	0.009	-5.589	0.000
rad	0.016	0.003	5.700	0.000
tax	-0.001	0.000	-4.322	0.000
ptratio	-0.036	0.006	-5.941	0.000
b	0.000	0.000	3.447	0.001
lstat	-0.030	0.002	-13.213	0.000

Notad que la edad sigue sin ser significativa. Pero esto puede deberse a que la relación es no lineal en edad.

Modelos de regresión lineal regularizada

Aquí vamos a considerar solamente un variante de regresión regularizada, LASSO, pero debéis recordar que también existe “ridge regression”. De hecho, hemos discutido que en determinadas situaciones puede ser preferible a LASSO.

A partir del análisis exploratorio de datos deberíamos decidir tanto qué variables considerare en el modelo como la forma de incluirla (no linealidad, interacciones, etc.). En el caso de LASSO, no es muy crucial puesto que el propio algoritmo se encargará de no seleccionar aquellas no relevantes, así que podemos incluir más variables de las que a priori consideremos más útiles. En otros algoritmos, esto puede tener efectos negativos sobre el funcionamiento del modelo (ej. “overfitting”).

En cualquier caso, existen varias cuestiones a tener en cuenta. En primer lugar, tenemos indicios de que puede ser conveniente transformar la variables dependiente en logaritmos. Esto implica que consideraremos estimar los modelos tanto para la variable en niveles como en logaritmos. También tenemos que considerar distintas variantes de transformaciones para las variables explicativas: ¿qué interacciones considerar? ¿qué formas de no linealidad: polinomios o discretización de variables discretas? Para variables categóricas, ¿reagrupamos categorías?

Por simplicidad, solo considero unas pocas especificaciones (diferentes combinaciones de variables); vosotros debéis considerar más, aunque quizás no tengáis que reportar todas. Lo importante es justificar por qué se ha decidido probar esas variantes (en particular, en función del análsis exploratorio) y por qué se han elegido las que consideréis. NOTAD que podemos partir directamente de un modelo LASSO muy general sin haber hecho ningún modelo de regresión lineal previo.

En particular, en los modelos LASSO solo mostraremos resultados para la variable dependiente en logaritmos; pero vosotros debéis probar y mostrar alguno con la variable dependiente en niveles. Sí consideraré algunas especificaciones usando polinomios o discretizando algunas variables continuas para capturar no linealidades y usando interacciones.

Además de especificar las variables dependiente y explicativas del modelo, la receta debe incluir un preprocesado de las variables adecuado a este algoritmo. No se pueden pasar factores a LASSO, por lo que debemos convertirlos en dummies con step_dummy, y las variables continuas se deben estandarizar.

recetaLASSO1 <- training(bostonPart) %>%            
  recipe(lmedv ~ chas + crim + zn + indus + nox + rm + 
           age + dis + rad + tax + ptratio + b + lstat) %>% 
  step_cut(b, breaks = c(0, 350, 395, 500), 
           include_outside_range = TRUE) %>% 
  step_normalize(all_predictors(), -all_nominal()) %>%
  step_poly(lstat, age, degree = 8) %>% 
  step_dummy(b, chas) %>% 
  step_interact(terms = ~ rm:(nox + starts_with("age_") + starts_with("lstat_"))) %>%  
  step_interact(terms = ~ nox:(crim + starts_with("age_") + starts_with("lstat_"))) %>% 
  step_interact(terms = ~ starts_with("b_"):(crim + starts_with("age_")))

recetaLASSO2 <- training(bostonPart) %>%            
                recipe(lmedv ~ chas + crim + zn + indus + nox + rm + 
                    age + dis + rad + tax + ptratio + b + lstat) %>% 
                step_cut(age, breaks = seq(0,100,10), 
                         include_outside_range = TRUE) %>% 
                step_cut(dis, breaks = c(0, 3, 10), 
                         include_outside_range = TRUE) %>% 
                step_cut(b, breaks = c(0, 350, 395, 500), 
                         include_outside_range = TRUE) %>% 
                step_normalize(all_predictors(), -all_nominal()) %>%
                step_poly(lstat, degree = 8) %>% 
                step_dummy(age, dis, b, chas)  %>% 
                step_interact(terms = ~ rm:(nox + starts_with("age_") + starts_with("lstat_"))) %>%  
                step_interact(terms = ~ nox:(crim + starts_with("age_") + starts_with("lstat_"))) %>% 
                step_interact(terms = ~ starts_with("b_"):(crim + starts_with("age_")))

Definimos los modelos, con el hiperparámetro de penalización para ajustar, y combinamos la receta y el modelo en un flujo.

modeloLASSO <- linear_reg(mode= "regression", penalty = tune()) %>%
                    set_engine("glmnet")

flujoLASSO1 <- workflow() %>%
  add_recipe(recetaLASSO1) %>%
  add_model(modeloLASSO)

flujoLASSO2<- workflow() %>%
  add_recipe(recetaLASSO1) %>%
  add_model(modeloLASSO)

La estimación del modelo implica ajustar el hiperparámetro $\lambda$ mediante validación cruzada en la muestra de entrenamiento. Necesitamos probar varios varios rangos de valores para el hiperparámetro.

set.seed(9753)
Boston_cv <- training(bostonPart) %>% vfold_cv(v=10)

LASSO_grid <- grid_regular(penalty(range = c(0, 1), trans = NULL),   
                          levels = 51)                     

set.seed(1)
lasso1_flujo_tuned <- flujoLASSO1 %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = LASSO_grid                          
                          ) 
lasso1_flujo_tuned %>% autoplot()

####
LASSO_grid <- grid_regular(penalty(range = c(0, 0.02), trans = NULL),   
                          levels = 51)                     

set.seed(1)
lasso1_flujo_tuned <- flujoLASSO1 %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = LASSO_grid                          
                          ) 

lasso1_flujo_tuned %>% autoplot()

Figura 11. Ajuste de lamdba, modelo 1

lasso1_flujo_tuned %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 11. Mejores lambdas en el Modelo 1 de LASSO") %>% kable_paper("hover")

Tabla 11. Mejores lambdas en el Modelo 1 de LASSO
penalty	.metric	.estimator	mean	n	std_err	.config
0.0036	rmse	standard	0.1689	10	0.0123	Preprocessor1_Model10
0.0040	rmse	standard	0.1690	10	0.0119	Preprocessor1_Model11
0.0044	rmse	standard	0.1692	10	0.0116	Preprocessor1_Model12
0.0032	rmse	standard	0.1692	10	0.0129	Preprocessor1_Model09
0.0048	rmse	standard	0.1695	10	0.0114	Preprocessor1_Model13

# lambda1 <- 0 #
lambda1 <- lasso1_flujo_tuned %>% select_best(metric = "rmse")

Notad que muestro dos búsquedas de valores, refinando en la segunda la zona que aparecía como más probable en la primera. En general, se deberían probar varios rangos para el hiper-parámetro buscando una forma de U para el valor con mínimo error; se puede empezar con rangos amplios y luego buscar valores en un rango más fino, entorno al valor donde se ve el mínimo. NO está claro que queráis mostrarlos todos o incluso que queráis mostrar más de uno.

LASSO_grid <- grid_regular(penalty(range = c(0, 0.2), trans = NULL),   
                          levels = 51)                     

set.seed(1)
lasso2_flujo_tuned <- flujoLASSO2 %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = LASSO_grid                          
                          ) 

lasso2_flujo_tuned %>% autoplot()

####
LASSO_grid <- grid_regular(penalty(range = c(0, 0.02), trans = NULL),   
                          levels = 51)                     

set.seed(1)
lasso2_flujo_tuned <- flujoLASSO2 %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = LASSO_grid                          
                          ) 

lasso2_flujo_tuned %>% autoplot()

lasso2_flujo_tuned  %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 12. Mejores lambdas en el Modelo 2 de LASSO") %>% kable_paper("hover")

Tabla 12. Mejores lambdas en el Modelo 2 de LASSO
penalty	.metric	.estimator	mean	n	std_err	.config
0.0036	rmse	standard	0.1689	10	0.0123	Preprocessor1_Model10
0.0040	rmse	standard	0.1690	10	0.0119	Preprocessor1_Model11
0.0044	rmse	standard	0.1692	10	0.0116	Preprocessor1_Model12
0.0032	rmse	standard	0.1692	10	0.0129	Preprocessor1_Model09
0.0048	rmse	standard	0.1695	10	0.0114	Preprocessor1_Model13

# lambda2 <- 0 
lambda2 <- lasso2_flujo_tuned %>% select_best(metric = "rmse")

Vamos a finalizar los modelos y ver sus métricas en la muestra de prueba.

LASSO_final <- list() 
LASSO_final[[1]] <- flujoLASSO1 %>% 
                    finalize_workflow(lambda1) %>% 
                    last_fit(bostonPart)  %>% 
                    collect_metrics()

LASSO_final[[2]] <- flujoLASSO2 %>% 
                    finalize_workflow(lambda2) %>% 
                    last_fit(bostonPart)  %>% 
                    collect_metrics()

LASSO_final %>% bind_rows(.id = "modelo") %>% 
  pivot_wider(names_from = .metric, values_from=.estimate) %>% 
  select(-.estimator) %>% 
  kbl(digits = 4, caption = "Tabla 13. Métricas de LASSO") %>% kable_paper("hover")

Tabla 13. Métricas de LASSO
modelo	.config	rmse	rsq
1	Preprocessor1_Model1	0.2411	0.67
2	Preprocessor1_Model1	0.2411	0.67

En este caso, ambos modelos tienen un comportamiento predictivo muy similar. Esto puede deberse a que ambos son muy similares, salvo por usar polinomios en lugar de discretizaciones. No es raro que ambos tengan un comportamiento predictivo similar, puesto que ambas formas aproximan (aunque de forma ligeramente distinta) relaciones no lineales.

kNN

Consideremos modelos de k vecinos. Este método es suficientemente flexible para considerar posibles no linealidades en todas las variables. Vamos a ajustar el número de vecinos por validación cruzada con el precio y su logaritmo como variables dependientes. Podríamos considerar otras variantes como usar el absoluto para la distancia o usar más o menos variables explicativas en el modelo; también usar discretizaciones de variables continuas (aunque esto no es muy habitual ni útil en kNN).

NOTA: como vemos aquí, conviene empezar con un rango amplio de valores para el hiperparámetro, sin probar cada valor del rango (ej., entre 1 y 21 aumentando valores de 2 en 2), y luego refinar en la zona donde parece estar el mínimo.

knn1_receta <- training(bostonPart) %>%            
  recipe(medv ~ chas + crim + zn + indus + nox + rm + 
           age + dis + rad + tax + ptratio + b + lstat) %>% 
  step_normalize(all_predictors(), -all_nominal()) %>%
  step_dummy(all_nominal_predictors())
  
knn1_modelo <- nearest_neighbor(mode = "regression",
                  neighbors = tune(), dist_power = 2) %>% 
                set_engine("kknn")

knn1_flujo <- workflow() %>%
  add_recipe(knn1_receta) %>%
  add_model(knn1_modelo)

knn_grid <- grid_regular(neighbors(range = c(1, 21), trans = NULL),   
                          levels = 11)                     

set.seed(1)
knn1_flujo_tuned <- knn1_flujo %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = knn_grid                          
                          ) 

knn1_flujo_tuned  %>% autoplot()

Figura 13. Ajuste de k vecinos, modelo 1

knn_grid <- grid_regular(neighbors(range = c(1, 11), trans = NULL),   
                          levels = 11)                     
###########
set.seed(1)
knn1_flujo_tuned <- knn1_flujo %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = knn_grid                          
                          ) 

knn1_flujo_tuned  %>% autoplot()

knn1_flujo_tuned  %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 14. k-vecinos óptimo en el Modelo 1") %>% kable_paper("hover")

Tabla 14. k-vecinos óptimo en el Modelo 1
neighbors	.metric	.estimator	mean	n	std_err	.config
4	rmse	standard	4.0733	10	0.3965	Preprocessor1_Model04
3	rmse	standard	4.0959	10	0.4083	Preprocessor1_Model03
5	rmse	standard	4.1115	10	0.3893	Preprocessor1_Model05
6	rmse	standard	4.1555	10	0.3876	Preprocessor1_Model06
7	rmse	standard	4.1946	10	0.3897	Preprocessor1_Model07

k1 <- knn1_flujo_tuned %>% select_best(metric = "rmse")

knn2_receta <- training(bostonPart) %>%            
  recipe(lmedv ~ chas + crim + zn + indus + nox + rm + 
           age + dis + rad + tax + ptratio + b + lstat) %>% 
  step_normalize(all_predictors(), -all_nominal()) %>%
  step_dummy(all_nominal_predictors())

knn2_flujo <- workflow() %>%
  add_recipe(knn2_receta) %>%
  add_model(knn1_modelo)

knn_grid <- grid_regular(neighbors(range = c(1, 11), trans = NULL),   
                          levels = 11)  

set.seed(1)
knn2_flujo_tuned <- knn2_flujo %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = knn_grid                          
                          ) 

knn2_flujo_tuned  %>% autoplot()

Figura 14. Ajuste de k vecinos, modelo 2

knn2_flujo_tuned  %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 15. k-vecinos óptimo en el Modelo 2") %>% kable_paper("hover")

Tabla 15. k-vecinos óptimo en el Modelo 2
neighbors	.metric	.estimator	mean	n	std_err	.config
5	rmse	standard	0.1722	10	0.0099	Preprocessor1_Model05
6	rmse	standard	0.1726	10	0.0104	Preprocessor1_Model06
4	rmse	standard	0.1726	10	0.0096	Preprocessor1_Model04
7	rmse	standard	0.1733	10	0.0107	Preprocessor1_Model07
8	rmse	standard	0.1745	10	0.0110	Preprocessor1_Model08

k2 <- knn2_flujo_tuned %>% select_best(metric = "rmse")

El número óptimo de vecinos es

neighbors	.config
4	Preprocessor1_Model04

neighbors	.config
5	Preprocessor1_Model05

knn_final <- list() 
knn_final[[1]] <- knn1_flujo %>% 
                    finalize_workflow(k1) %>% 
                    last_fit(bostonPart)  %>% 
                    collect_metrics()

knn_final[[2]] <- knn2_flujo %>% 
                    finalize_workflow(k2) %>% 
                    last_fit(bostonPart)  %>% 
                    collect_metrics()

knn_final %>% bind_rows(.id = "modelo") %>% 
  pivot_wider(names_from = .metric, values_from=.estimate) %>% 
  select(-.estimator) %>% 
  kbl(digits = 4, caption = "Tabla 16. Métricas de kNN") %>% kable_paper("hover")

Tabla 16. Métricas de kNN
modelo	.config	rmse	rsq
1	Preprocessor1_Model1	4.2681	0.7289
2	Preprocessor1_Model1	0.2031	0.7350

Comparando estos resultados con los de la Tabla para LASSO, queda claro que kNN tiene una mejor capacidad predictiva, tanto para la variable de precio como usando el logaritmo del precio. Sin embargo, estos modelos no nos dicen mucho sobre qué factores afectan a esa predicción (no tenemos ni modelo paramétrico y contrastes, ni representación gráficas o medidas de importancia).

Árboles de regresión

Pasamos a considerar modelos de árboles de decisión. Notad que para ajustar el coste de complejidad puede ser necesario varios rangos. Dado que hemos visto que sistemáticamente el modelo con el logaritmo del precio predice mejor, nos centramos en ese modelo (aunque esto es obviamene una elección y quizás deberían probarse más).

arbol2_receta <- training(bostonPart) %>%            
  recipe(lmedv ~ chas + crim + zn + indus + nox + rm + 
           age + dis + rad + tax + ptratio + b + lstat) 

arbol2_modelo <-  decision_tree(mode = "regression", 
                                     cost_complexity = tune()) %>% 
                          set_engine("rpart") 

arbol2_flujo <- workflow() %>%
  add_recipe(arbol2_receta) %>%
  add_model(arbol2_modelo)

arbol_grid <- grid_regular(cost_complexity(range = c(0, 0.01), trans = NULL),   
                          levels = 11)                     
set.seed(1)
arbol2_flujo_tuned <- arbol2_flujo %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = arbol_grid                          
                          ) 

arbol2_flujo_tuned  %>% autoplot()

arbol2_flujo_tuned  %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 17. Coste de complejidad  en el Modelo 1") %>% kable_paper("hover")

Tabla 17. Coste de complejidad en el Modelo 1
cost_complexity	.metric	.estimator	mean	n	std_err	.config
0.002	rmse	standard	0.1898	10	0.0119	Preprocessor1_Model03
0.001	rmse	standard	0.1898	10	0.0122	Preprocessor1_Model02
0.000	rmse	standard	0.1899	10	0.0123	Preprocessor1_Model01
0.003	rmse	standard	0.1935	10	0.0123	Preprocessor1_Model04
0.004	rmse	standard	0.1971	10	0.0117	Preprocessor1_Model05

cost2 <- arbol2_flujo_tuned %>% select_best(metric = "rmse")

arbol_final <- arbol2_flujo %>% 
                    finalize_workflow(cost2) %>% 
                    last_fit(bostonPart)  

arbol_final %>%  collect_metrics() %>%  
  pivot_wider(names_from = .metric, values_from=.estimate) %>% 
  select(-.estimator) %>% 
  kbl(digits = 4, caption = "Tabla 18. Métricas del árbol") %>% kable_paper("hover")

Tabla 18. Métricas del árbol
.config	rmse	rsq
Preprocessor1_Model1	0.2064	0.7394

Vemos que el mejor árbol de regresión es mejor que LASSO y similar a kNN (en logaritmos). Además podemos ver el flujo del árbol para hacernos una idea de qué factores están afectando a la predicción

arb <- arbol_final$.workflow[[1]] %>%   
  extract_fit_parsnip() 

rpart.plot(arb$fit)

Figura 16a. Mejor árbol de regresión, modelo 2

También podemos recurrir a medidas de importancia de las variables.

arbol_final$.workflow[[1]] %>%   
  extract_fit_parsnip() %>% 
  vip(num_features = 14)

Figura 16b. Importancia en el mejor árbol de regresión, modelo 2

Estos resultado son interesantes, más allá de su capacidad predictiva, porque son informativos son las interacciones entre distintas características. Como sabemos y podemos ver, los árboles trabajan discretizando las variables continuas discretizadas. Pero además implican distintas interacciones entre rangos de las variables para el resultado final. Por ejemplo, vemos que en los niveles superiores del árbol tenemos distintos niveles de población desfavorecida en el barrio, lstat; de hecho, vemos también (como era esperable) que es la variable con mayor valor de la importancia. Pero según los distintos valores de lstat importan diferentes características: cuando lstat es bajo ($<10$) el precio viene determinado por el número de habitaciones de la vivienda, mientras que cuando es alto por la cantidad de contaminación o la proporición de gente de color en el vecindario.

Random Forests

Intentamos ver si un modelo de “random forests” puede mejorar al árbol de regresión u otros modelos.

rf2_receta <- arbol2_receta 

rf2_modelo <-  rand_forest(mode = "regression",
                                       mtry = tune(), trees = 100) %>% 
                          set_engine("ranger", importance = "impurity")

rf2_flujo <- workflow() %>%
  add_recipe(rf2_receta) %>%
  add_model(rf2_modelo)

rf_grid <- grid_regular(mtry(range = c(1, 12), trans = NULL),   
                          levels = 12)                     

rf2_flujo_tuned <- rf2_flujo %>% 
                        tune_grid(
                          resamples = Boston_cv,
                          metrics   = metric_set(rmse),
                          grid      = rf_grid                          
                          ) 

rf2_flujo_tuned %>% autoplot()

rf2_flujo_tuned  %>% show_best(metric = "rmse") %>% 
  kbl(digits = 4, caption = "Tabla 19. Ajuste de Random Forest  en el Modelo 2") %>% kable_paper("hover")

Tabla 19. Ajuste de Random Forest en el Modelo 2
mtry	.metric	.estimator	mean	n	std_err	.config
7	rmse	standard	0.1363	10	0.0086	Preprocessor1_Model07
9	rmse	standard	0.1368	10	0.0085	Preprocessor1_Model09
6	rmse	standard	0.1378	10	0.0085	Preprocessor1_Model06
8	rmse	standard	0.1383	10	0.0095	Preprocessor1_Model08
12	rmse	standard	0.1391	10	0.0091	Preprocessor1_Model12

mtry2 <- rf2_flujo_tuned %>% select_best(metric = "rmse")

El número óptimo de variables a usar (aleatoriamente) en cada nodo son

mtry	.config
7	Preprocessor1_Model07

Por último finalizamos el modelo

rf_final <- rf2_flujo %>% 
                    finalize_workflow(mtry2) %>% 
                    last_fit(bostonPart)  

rf_final %>%  collect_metrics() %>%  
  pivot_wider(names_from = .metric, values_from=.estimate) %>% 
  select(-.estimator) %>% 
  kbl(digits = 4, caption = "Tabla 20. Métricas de RF") %>% kable_paper("hover")

Tabla 20. Métricas de RF
.config	rmse	rsq
Preprocessor1_Model1	0.1619	0.824

El modelo de “Random Forest” incluso mejora claramente a “kNN”. En este caso, podemos ver qué variables han resultados mejores predictores del precio.

rf_final$.workflow[[1]] %>%   
  extract_fit_parsnip() %>% 
  vip(num_features = 14)

Figura 18. Importancia en el mejor RF, modelo 2

Este gráfico nos indica que la variable más importante (con diferencia) es la relacionada con la “pobreza” de la zona en la que se encuentra la casa. Después nos encontramos una característica de la casa en sí (número de habitaciones, aunque realmente es la media de la zona) y la criminalidad de la zona. No existe un criterio formal para decidir en que punto parar determinar qué variables son importantes. Quizás en este caso se puede “argumentar” que hay una caída mayor de la importancia después de la contaminación de la zona, la calidad de las escuelas y la distancia al centro.

Conclusión

Uno de los objetivos de este informe era determinar los atributos del vecindario que mejor explicaban (no necesariamente causan) la variación en los precios de las viviendas. Los mejores modelos son regresión lineal, que incorpora relaciones no lineales, y “random forest” que también establece una relación altamente no lineal entre el precio y los factores observados. También ofrecen herramientas para interpretar los resultados. Ambos modelos coinciden en la importancia de que el barrio esté en una zona “deprimida” (mucho porcentaje de población con bajo estatus social) en la determinación del precio. Con mayor o menor importancia, también coinciden en la relevancia de la criminalidad en la zona, el número medio de habitaciones y el ratio de alumnos por profesor. Algunas características como la edad o una localización cercana al río no parecen ser importantes; esto es una diferencia respecto al análisis de covariación, donde sí se apreciaba cierta correlación. Esto nos indica que esas características están relacionadas con otras que sí son más relevantes en la fijación de precios. Los resultados del modelo de árboles también son relevantes para aprender sobre las interacciones entre características para determinar el precio de las casas. Debéis usar estos resultados para reestimar modelos más interpretables (como regresión lineal o árboles) para hacer comentarios específicos sobre el impacto de cambios concretos de variables particularmente relevantes sobre el precio de las casas. PENSAD en qué cuestiones pueden ser interesantes si estáis trabajando en una inmobiliaria o en una empresa de tasación de viviendas, en este caso, o en una empresa del sector apropiado para los datos que elijáis.

Por otro lado estos resultados también se podrían utilizar para predecir el precio de una nueva casa (o saber si una en el mercado está sobrevalorada).

rf_final$.workflow[[1]] %>% 
  predict(new_data = tibble(chas = "Yes", crim = 0.03,  zn = 0,
                            indus = 7, nox = 0.5, rm = 2,
                            age = 10,  dis = 6, rad = 5,  
                            tax = 290, ptratio = 20, b = 390, lstat=9)) %>% 
  exp()

.pred
20.58078

Notad que debemos deshacer el logaritmo para dar una predicción en unas unidades con sentido (miles de dólares). Recordad que también podemos ofrecer un intervalo de confianza para la predicción.

Fijaos que en otros contextos o con otros datos podéis tener un caso más obvio que predecir (o no). Y podéis usar intervalos de confianza en lugar de predicción puntual para determinar si una casa está sobrevalorada o infravalorada. En cualquier caso, debéis elegir y justificar los valores elegidos para predecir: en este caso, una casa observada en los datos con particular interés, una casa común en la zona. En el proyecto que elijáis, también deben existir razones para que ese caso sea relevante para predecir.

Recordad que también podemos estar interesados en usar el modelo para interpetar el efecto de alguna o algunas variables. Nuevamente, tendréis que justificar por qué esas variables son interesantes y cómo se interpretan los resultados. Como hemos visto, podemos considerar la información que LASSO y/o “Random Forest” nos dan sobre las variables para estimar modelos más interpretables y extraer conclusiones.

Apéndice A

Como se ha comentado, nuestro conjunto de datos dispone de 506 observaciones y 14 variables. A continuación se presenta una breve descripción de cada variable:

Variable	Descripción
CRIM	tasa de criminalidad per cápita por ciudad
ZN	proporción de terreno residencial dividido en
	zonas para lotes de más de 25,000 pies cuadrados.
INDUS	proporción de acres de negocios no comerciales por zona
CHAS	variable “dummy”del río Charles
	(1 si la zona limita con el río; si no, 0)
NOX	concentración de óxidos nítricos (partes por 10 millones)
RM	número medio de habitaciones por vivienda
AGE	proporción de casas ocupadas por sus propietarios
	construidas antes de 1940.
DIS	distancia media ponderada a cinco centros de empleo de Boston
RAD	índice de accesibilidad a carreteras radiales
TAX	impuesto sobre el valor total de la propiedad
	(por cada $10,000)
PTRATIO	número de alumnos por docente en los colegios de la zona (“town”)
B	1000(Bk - 0.63)^2 donde Bk es la proporción de población de color
	en la zona (0.63 es la media en la ciudad)
LSTAT	porcentaje de población desfavorecida (bajo estatus social)
MEDV	valor mediano de las viviendas ocupadas por sus propietarios
	(en miles de dólares)

Notas

Todo esta discusión sobre modelos de regresión lineal aplica a regresión logística para problemas de clasificación.↩︎