Prediciendo precios de casas de Ames: análisis exploratorio de datos

Configuración

%%capture
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import warnings
from scipy.stats import norm, probplot
from sklearn.datasets import fetch_openml
from utils import plot_univariate, plot_bivariate
from IPython.display import display
pd.set_option('display.max_columns', 200)
warnings.filterwarnings('ignore')
sns.distributions._has_statsmodels=False # Estimating density with scipy

El conjunto de datos contiene 79 valores por cada observación, 80 en el conjunto de entrenamiento considerando la variable objetivo. Hay 1460 observaciones en el conjunto de entrenamiento y una menos en el conjunto que se utiliza para evaluar el análisis, es decir, predecir la variable objetivo, que es SalePrice. La única columna que no se utilizará es Id.

X_train = fetch_openml(name="house_prices", as_frame=True,
                       data_home='data')['frame']
X_train.drop(['Id'], axis=1, inplace=True)
for col in ['MSSubClass', 'OverallQual', 'OverallCond']:
    X_train[col] = X_train[col].astype('int')
display(X_train.head(2), X_train.shape)

(1460, 80)

Las variables de entrada se han clasificado por sus niveles de medición (de razón, de intervalo, ordinal, nominal):

verbal = [f for f in X_train.columns if X_train.dtypes[f] == 'object']
numerical = [f for f in X_train.columns if X_train.dtypes[f] != 'object']
qualitative = verbal + ['MSSubClass', 'OverallQual', 'OverallCond']
quantitative = [
    n for n in numerical if n not in qualitative and n != 'SalePrice']
interval = ['GarageYrBlt', 'YearBuilt', 'YearRemodAdd', 'YrSold', 'MoSold']
ratio = [c for c in quantitative if c not in interval and c != 'SalePrice']
ratio.sort()
ordinal = ['LotShape', 'Utilities', 'LandSlope', 'ExterQual', 'OverallQual', 'OverallCond',
           'ExterCond', 'BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2',
           'HeatingQC', 'Electrical', 'KitchenQual', 'Functional', 'FireplaceQu',
           'GarageFinish', 'GarageQual', 'GarageCond', 'PavedDrive', 'PoolQC', 'Fence']
ordinal.sort()
nominal = [n for n in qualitative if n not in ordinal]
nominal.sort()

La respuesta: SalePrice

Los gráficos del lado izquierdo de la siguiente figura muestran la distribución de todos los precios que utilizaré para supervisar el análisis y la gráfica de probabilidad de todas estas muestras.

Podemos ver asimetría positiva y curtosis, hay una cola larga a la derecha y la distribución es muy puntiaguda alrededor de la media, por lo que nuestra respuesta no sigue una distribución normal (mostrada en negro).

Sin embargo, esto se rectifica parcialmente después de una transformación logarítmica, como se puede ver en las gráficas de probabilidad y distribución del lado derecho: la distribución sigue de cerca la diagonal en la gráfica de probabilidad correspondiente y la densidad está más cerca de la forma normal.

sns.set(font_scale=1.35, style='dark')
fig = plt.figure(figsize=(17, 9))
plt.subplots_adjust(hspace=.3)
title = ['SalePrice', 'log(SalePrice)']
for i, y in enumerate([X_train['SalePrice'], np.log(X_train['SalePrice'])]*2, start=1):
    ax = plt.subplot(2, 2, i)
    if i < 3:
        sns.distplot(y, fit=norm, ax=ax)
        ax.legend(['Normal distribution', 'Skewness: {:.3}\nKurtosis: {:.3}'.format(y.skew(),
                                                                                    y.kurt())],
                  loc='best')
        ax.set_title(title[i-1], fontdict={'fontweight': 'bold'})
        ax.set_xlabel(title[i-1])
    else:
        probplot(y, plot=ax)
        ax.get_lines()[0].set_markerfacecolor('#7b9ece')
plt.show()

X_train[['SalePrice']].describe()

Los predictores

Primero se realizará un análisis gráfico, mirando los predictores uno por uno, y luego las gráficas de los predictores como una función de la respuesta. Las cifras se clasifican por nivel de medida: de razón, de intervalo, ordinal y nominal. Utilizarés estas dos funciones para trazar los gráficos.

Análisis univariante

Este análisis es útil para descubrir algunas particularidades en los datos:

• La gran mayoría de las casas tienen menos de 2000 pies cuadrados de superficie habitable y menos de 20000 pies cuadrados de área de lote.
• La parcela de distribución de TotalBsmtSF muestra que la mayoría de las casas tienen sótano.
• La mayoría de las casas tienen entre 5 y 7 habitaciones sobre rasante, sin embargo, hay casas con hasta 14 habitaciones sobre rasante.
• Hay algunas casas con más de una cocina sobre rasante y algunas sin dormitorios.
• La mayoría de las propiedades tienen 1 o 2 baños completos, pero hay algunas con 3 y otras sin baños completos.
• La mayoría de las propiedades tienen suficiente espacio para hasta dos automóviles en el garaje.
• La mayoría de las casas no tienen piscina ni porche.

# Ratio-measured variables
X_train[ratio].describe()

# add '3SsnPorch', 'OpenPorchSF', 'ScreenPorch' below to check the rest of the variables with values for Porch
ratio_vars = ['BedroomAbvGr', 'EnclosedPorch', 'KitchenAbvGr', 'FullBath', 'GarageCars',
              'GrLivArea', 'LotArea', 'TotalBsmtSF', 'TotRmsAbvGrd']
plot_univariate(X_train, ratio_vars, 2.5, 3, 5, 2, sns.distplot)

plot_univariate(X_train, ratio_vars+['PoolArea'], 2, 3, 4, 2, sns.boxplot)

• Hay observaciones en YearRemodAdd desde 1950 en adelante y, dada la enorme diferencia en cantidad de este año a los siguientes cincuenta años, es probable que haya un límite inferior en sus valores. Además, como las propiedades nuevas probablemente no necesiten remodelación y YearRemodAdd no contiene valores que faltan, muchos de los valores, aparte de los de las casas con YearRemodAdd de 1950, pueden haber sido imputados de alguna manera prescriptiva.

# Interval-measured variables
X_train[interval].describe()

interval_vars = ['YearBuilt', 'YearRemodAdd', 'GarageYrBlt']
plot_univariate(X_train, interval_vars, 1.8, 3, 3, 2, sns.distplot, bins=50)

• Obviamente, existe una estrecha relación entre predictores como GarageYrBlt y YearBuilt, o GarageArea y GarageCars, sin embargo, dado que la mayoría de las casas tienen garaje y la distribución de GarageYrBlt y YearBuilt es ligeramente diferente, es probable que el garaje en muchas propiedades se haya construido más tarde o que la unidad de garaje y la casa sean propiedades independientes y, en consecuencia, lo más probable es que ambas se hayan construido en años diferentes.

print("Number of houses with garage:",
      X_train[X_train['GarageArea'] > 0].shape[0],
      "\nNumber of houses built the same year as the garage:",
      X_train[(X_train['GarageArea'] > 0) &
              (X_train['YearBuilt'] == X_train['GarageYrBlt'])].shape[0])

Number of houses with garage: 1379 
Number of houses built the same year as the garage: 1089

• Se han vendido muchas más casas durante la temporada de verano. Además, el número de casas vendidas de un año a otro desde 2006 no ha cambiado mucho excepto en 2010, sin embargo, esto se debe a que en el conjunto de datos hay registros de casas que se han vendido solo hasta julio de 2010, como podemos ver a continuación:

interval_vars2 = ['MoSold', 'YrSold']
plot_univariate(X_train, interval_vars2, 1.4, 2, 3, 2, sns.barplot)

display(pd.DataFrame(X_train['YrSold'].value_counts()))
print("Last month with records in 2010: ",
      X_train[X_train['YrSold'] == 2010]['MoSold'].max())

Last month with records in 2010:  7

# Nominal-measured variables
nominal_vars = ['BldgType', 'CentralAir', 'Heating', 'HouseStyle', 'LandContour', 'LotConfig',
                'MSZoning', 'MSSubClass', 'Neighborhood', 'SaleCondition', 'SaleType']
plot_univariate(X_train, nominal_vars, 1.6, 3, 3, 2,
                sns.barplot, labels_thresh=9, rotation=60)

• La gran mayoría de viviendas tienen aire acondicionado y también generadores de aire caliente de tiro forzado a gas como tipo de calefacción.
• La mayoría de las viviendas se han construido en una zona residencial de baja densidad.

# Ordinal-measured variables
X_train['OverallCond'] = X_train['OverallCond'].astype('category')
X_train['OverallQual'] = X_train['OverallQual'].astype('category')
frame = X_train[ordinal].describe()
frame.index = ["count", 'unique', 'mode', 'freq']
display(frame)
X_train['OverallCond'] = X_train['OverallCond'].astype('int')
X_train['OverallQual'] = X_train['OverallQual'].astype('int')

ordinal_vars = ['BsmtCond', 'BsmtQual', 'ExterQual', 'ExterCond',  'FireplaceQu', 'GarageCond',
                'GarageQual', 'HeatingQC', 'KitchenQual', 'LandSlope', 'OverallCond',
                'OverallQual']
frame_ordinal = X_train[ordinal_vars].copy()
Bsmt_Fireplace_Garage_mapping = {
    'Ex': 5, 'Gd': 4, 'TA': 3, 'Fa': 2, 'Po': 1, np.nan: 0}
Exter_Heating_Kitchen_mapping = {'Ex': 5, 'Gd': 4, 'TA': 3, 'Fa': 2, 'Po': 1}
LandSlope_mapping = {'Gtl': 2, 'Mod': 1, 'Sev': 0}
frame_ordinal['BsmtCond'] = frame_ordinal['BsmtCond'].map(
    Bsmt_Fireplace_Garage_mapping)
frame_ordinal['BsmtQual'] = frame_ordinal['BsmtQual'].map(
    Bsmt_Fireplace_Garage_mapping)
frame_ordinal['ExterQual'] = frame_ordinal['ExterQual'].map(
    Exter_Heating_Kitchen_mapping)
frame_ordinal['ExterCond'] = frame_ordinal['ExterCond'].map(
    Exter_Heating_Kitchen_mapping)
frame_ordinal['FireplaceQu'] = frame_ordinal['FireplaceQu'].map(
    Bsmt_Fireplace_Garage_mapping)
frame_ordinal['GarageCond'] = frame_ordinal['GarageCond'].map(
    Bsmt_Fireplace_Garage_mapping)
frame_ordinal['GarageQual'] = frame_ordinal['GarageQual'].map(
    Bsmt_Fireplace_Garage_mapping)
frame_ordinal['HeatingQC'] = frame_ordinal['HeatingQC'].map(
    Exter_Heating_Kitchen_mapping)
frame_ordinal['KitchenQual'] = frame_ordinal['KitchenQual'].map(
    Exter_Heating_Kitchen_mapping)
frame_ordinal['LandSlope'] = frame_ordinal['LandSlope'].map(LandSlope_mapping)
plot_univariate(frame_ordinal, ordinal_vars, 1.6, 3, 3, 2, sns.barplot)

• La mayoría de las propiedades son casi planas (pendiente suave del terreno). Obviamente, esto puede deberse al barrio o área donde se han construido. Sin embargo, puede haber algunas excepciones que pueden crear distinciones entre unas casas y otras en la misma área.

Algunas otras posibles distinciones son las casas que tienen diferentes calificaciones en OverallQual o en OverallCond o diferente LotArea dentro del mismo" vecindario.

Dado que estas diferencias pueden resultar en un aumento o disminución del precio, es posible que necesitemos crear algunas variables para explicar esta información en caso de que nuestro modelo no pueda capturar o explicar esta proporción de varianza en la respuesta.

Por ejemplo, valores atípicos y patrones en las gráficas de los residuos una vez que se ajusta el modelo, ya que esto podría indicar una deficiencia en el análisis, como la falta de una variable.

Análisis bivariante

Las gráficas siguientes nos ayudan a identificar algunas observaciones que están fuera de la tendencia general del resto de los datos en cada uno de los espacios de predictores y objetivo.

Por ejemplo, sabemos por el diagrama de caja de GrLivArea que la misma contiene observaciones sobre el bigote superior, es decir, con un valor extremo en el espacio del predictor, pero ahora, como podemos ver a continuación, podemos afirmar que dos de de estos cuatro supuestos valores atípicos son verdaderamente valores atípicos, en el sentido de que no siguen la tendencia general de los datos.

Otro ejemplo claro es la propiedad con más de 6000 pies cuadrados de sótano en TotalBsmtSF, que por el precio podría ser uno de los valores atípicos previamente definidos.

Por otro lado, las otras dos observaciones en GrLivArea con más de 4000 pies cuadrados están en concordancia con el resto de puntos o, al menos, podríamos decir que siguen la tendencia del resto de casas.

# Ratio-measured variables and 'SalePrice'
plot_bivariate(X_train, 'SalePrice', ratio_vars, 2.8, 3, 5, 2, plt.scatter)

# Interval-measured variables and 'SalePrice'
plot_bivariate(X_train, 'SalePrice', interval_vars, 2.2, 2, 6, 2, plt.scatter)

Algunas casas tienen los mismos valores en YearRemodAdd como en YearBuilt por lo que aún no han sido remodeladas, es decir, algunos de sus valores podrían haber sido introducidos prescriptivamente. Esto puede deberse a que la propiedad se acaba de construir o que se encuentra en buenas condiciones.

Sin embargo, puede que no siempre sea el caso, como podemos ver en la tabla a continuación. Hay 39 propiedades en el conjunto de entrenamiento con un OverallCond inferior a 5 que se construyeron el mismo año en el que se especifica que fueron remodeladas:

pd.DataFrame(X_train[(X_train['YearRemodAdd'] == X_train['YearBuilt']) &
                     (X_train['OverallCond'] < 5)]['OverallCond'].value_counts())

La creación de algunas características podría capturar información que podría explicar la diferencia de tiempo entre YearBuilt y YrSold, o entre YearBuilt y YearRemodAdd, o entre YearRemodAdd y YrSold. Es decir, podría haber un patrón que podría ayudar a predecir el precio.

Como podemos ver a continuación tanto en el gráfico de la izquierda como en la tabla, el mes con la mediana más baja de SalePrice es abril. Enero, agosto y noviembre son los meses en los que hay más variabilidad en el precio de las casas en venta.

Aunque solo existen 1460 observaciones en este conjunto de datos, podemos decir que la primera mitad del año y julio es cuando hay más propiedades alejadas del grueso de los datos, i.e., quizás esto represente el tipo de casas que no se suelen vender durante el resto del año.

Además, octubre no parece el mejor mes para los agentes inmobiliarios. Tampoco lo fue 2008 en general.

Aunque la mediana de 2010 es incluso más baja que la mediana de 2008, solo hay casas en este conjunto hasta julio de 2010, como se mostró anteriormente.

plot_bivariate(X_train, 'SalePrice', interval_vars2, 1.6, 2, 4, 2,
               sns.boxplot)

pd.DataFrame(X_train.groupby(['MoSold'])
             ['SalePrice'].median().sort_values().head(3))

A continuación se muestran diagramas de caja de algunos predictores codificados ordinales frente a SalePrice . La codificación ayuda a comprender cómo influye el orden de los predictores en el precio de venta.

# Ordinal predictors and Sale Price
# Please refer to the mapping in univariate ordinal plot in order to check the x-labels
frame_ordinal = pd.concat([frame_ordinal, X_train['SalePrice']], axis=1)
plot_bivariate(frame_ordinal, 'SalePrice', ordinal_vars, 2.2, 3, 5, 2,
               sns.boxplot)

La mayoría de las veces existe una correlación positiva entre el precio y los predictores ordinales, sin embargo, si tomamos en cuenta solo los predictores que miden las condiciones, más no siempre es mejor, o una mejor condición puede no traducirse en un precio más alto.

De hecho, el 73,84% de las primeras 730 casas más caras del conjunto (la mitad del conjunto de datos) tienen un OverallCond de cinco de nueve:

series_cond = X_train.sort_values(by='SalePrice',
                                  ascending=False)['OverallCond'][:730].value_counts()
pd.DataFrame([series_cond, np.round(series_cond/series_cond.sum()*100, 2)],
             index=['n_houses', '%']).T

O también podemos decir que el 93% de las casas con un OverallQual de 9 o superior ha sido calificado con un OverallCond de 5 sobre 9:

series_cond = X_train[X_train['OverallQual'] > 8]['OverallCond'].value_counts()
pd.DataFrame([series_cond, np.round(series_cond/series_cond.sum()*100, 2)],
             index=['n_houses', '%']).T

Como podemos ver a continuación, y como se esperaba, el precio de las propiedades fluctúa mucho dependiendo del barrio y también en el tipo de vivienda (MSSubClass) entre otros. Por ejemplo, 1 pie cuadrado tiene una calificación muy diferente en StoneBr que en IDOTRR.

# Nominal predictors and Sale Price
plot_bivariate(X_train, 'SalePrice', nominal_vars, 2.6, 3, 5, 2, sns.boxplot,
               labels_thresh=7, rotation=60)  # rotation=70)

Podemos ver en la siguiente tabla los precios por pie cuadrado en promedio en los barrios más caros y más baratos (considerando solo el precio por pie cuadrado de la superficie habitable):

X_train['Living_Area_'] = X_train['GrLivArea'] + X_train['TotalBsmtSF']
dollar_sf = X_train.groupby(['Neighborhood'])[['Living_Area_',
                                               'SalePrice']].apply(lambda x: np.round(x.mean(), 2))
dollar_sf['$/SF'] = np.round(dollar_sf['SalePrice'] /
                             dollar_sf['Living_Area_'], 2)
display(dollar_sf.sort_values('$/SF', ascending=False).head(2))
display(dollar_sf.sort_values('$/SF', ascending=False).tail(2))

Si agrupamos por tipo de vivienda, así como de barrio, incluso podemos ver una diferencia mayor:

dollar_sf = X_train.groupby(['Neighborhood',
                             'MSSubClass'])[['Living_Area_',
                                             'SalePrice']].apply(lambda x: np.round(x.mean(), 2))
dollar_sf['$/SF'] = np.round(dollar_sf['SalePrice'] /
                             dollar_sf['Living_Area_'], 2)
display(dollar_sf.sort_values('$/SF', ascending=False).head(2))
display(dollar_sf.sort_values('$/SF', ascending=False).tail(2))

En promedio, el precio de un pie cuadrado de una vivienda de 2 pisos en Edwards es 5 veces más barato que en StoneBr. Esto es solo una aproximación a modo de ilustración, ya que muchas otras características del conjunto tienen una influencia, hasta cierto punto, en el precio del pie cuadrado.

Esto se explora más en este análisis supervisado.