1- Introducción
2- Caso de estudio y exploración de los datos
3- Procesamiento y limpieza de los datos
4- Ingeniería de características
5- Modelado, evaluación y predicción
Nuestra Comunidad
Votación: Temas de Interés en Ciencia de Datos con Python
In [1]:
# Librerías
import pandas as pd
import numpy as np
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV
from sklearn.metrics import mean_squared_error, make_scorer
from scipy import stats
from scipy.stats import norm, skew
from IPython.display import display
import matplotlib.pyplot as plt
import seaborn as sns
# Opciones
pd.set_option('display.float_format', lambda x: '%.3f' % x)
%matplotlib inline
In [2]:
# Para verificar los archivos contenidos en el directorio
from subprocess import check_output
print(check_output(["ls", "../python_regression/"]).decode("utf8"))
In [3]:
# Para leer los datos
train = pd.read_csv("train.csv")
print("train : " + str(train.shape))
test = pd.read_csv('test.csv')
print("test : " + str(test.shape))
In [4]:
## Despliega las primeras cinco filas del conjunto de datos de entrenamiento
train.head(5)
Out[4]:
In [5]:
## Describe los datos de entrenamiento
train.describe()
Out[5]:
In [6]:
# Verifica el número de observaciones y atributos
print("The train data size before dropping Id feature is : {} ".format(train.shape))
print("The test data size before dropping Id feature is : {} ".format(test.shape))
# Guarda la columna 'Id'
train_ID = train['Id']
test_ID = test['Id']
# Ahora prescinde de la columna 'Id', ya que es innecesaria para el proceso de predicción
train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)
# Verifica de nuevo el tamaño de los datos después de eliminar la variable 'Id'
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape))
print("The test data size after dropping Id feature is : {} ".format(test.shape))
La documentación de los datos de Ames Housing indica que existen valores extremos (outliers) en los datos de entrenamiento.
Vamos a explorar estos valores.
In [7]:
fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()
In [8]:
# Eliminamos los outliers
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)
# Verificamos de nuevo el gráfico
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()
Eliminar los valores extremos (outliers) no es siempre lo más apropiado. En nuestro caso procedimos a eliminar estos dos casos, debido a que eran valores muy elevados y que no hacían sentido (áreas extremadamente grandes a precios muy bajos).
Puede que existan otros valores extremos en los datos de entrenamiento. Sin embargo, removerlos todos puede afectar de mala manera nuestro ejercicio de modelado si también encontramos outliers en los datos de comprobación (test). Por eso, en lugar de removerlos todos, preferimos solo manejar aquellos casos que hagan nuestros modelos más robustos.
In [9]:
sns.distplot(train['SalePrice'] , fit=norm);
# Obtenemos los parámetros de ajuste utilizados por la función
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
# Graficamos la distribución empírica
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')
# QQ-plot para comprobar la normalidad
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()
In [10]:
# Empleamos la función log1p de numpy para obtener el log(1+x) de todos los elementos de la variable objetivo
train["SalePrice"] = np.log1p(train["SalePrice"])
# Verificamos la nueva distribución empírica
sns.distplot(train['SalePrice'] , fit=norm);
# Nuevos parámetros
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
# Nuevo gráfico después de la transformación
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')
# QQ-plot después de la transformación
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()
# Registramos la nueva variable target como y_train
y_train = train.SalePrice.values
In [11]:
ntrain = train.shape[0]
ntest = test.shape[0]
all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))
In [12]:
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(30)
Out[12]:
In [13]:
f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='0')
sns.barplot(x=all_data_na, y=all_data_na.index)
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)
Out[13]:
PoolQC : la descripción de datos dice que NA signifca "No Pool". Esto hace sentido, dado el enorme porcentaje de valores perdidos (+99%), y el hecho de que, en general, la mayoría de las viviendas no tienen piscina.
In [14]:
all_data["PoolQC"] = all_data["PoolQC"].fillna("None")
Casos similares.
In [15]:
all_data["MiscFeature"] = all_data["MiscFeature"].fillna("None")
all_data["Alley"] = all_data["Alley"].fillna("None")
all_data["Fence"] = all_data["Fence"].fillna("None")
all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")
for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
all_data[col] = all_data[col].fillna('None')
LotFrontage : Dado que el área de cada calle conectada a la casa, muy probablemente es similar para las otras casas del vecindario, podemos rellenar los valores perdidos con la mediana de LotFrontage del vecindario.
In [16]:
# Agrupamos por vecindario y completamos los valores perdidos con la mediana de LotFrontage para todos los vecindarios
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
lambda x: x.fillna(x.median()))
GarageYrBlt, GarageArea y GarageCars : Reemplazamos los datos perdidos con cero, ya que "no garage" equivale a que no hay carros en dicho garage.
In [17]:
for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
all_data[col] = all_data[col].fillna(0)
BsmtFinSF1, BsmtFinSF2, BsmtUnfSF, TotalBsmtSF, BsmtFullBath, BsmtHalfBath : los valores perdidos muy probablemente son cero por no tener basement.
In [18]:
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
all_data[col] = all_data[col].fillna(0)
BsmtQual, BsmtCond, BsmtExposure, BsmtFinType1, BsmtFinType2 : Para todos estos atributos categóricos relacionados con el basement, NaN significa que no hay basement.
In [19]:
for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
all_data[col] = all_data[col].fillna('None')
MasVnrArea y MasVnrType : NA muy probablemente indica que estas casas no tienen una chapa de albañilería. De manera que podemos reemplazar con cero el área y con None el tipo.
In [20]:
all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)
MSZoning (clasificación general de zonificación) : 'RL' es por mucho el valor más común. De manera que podemos reemplazar los valores perdidos con 'RL'.
In [21]:
all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])
Utilities : Para este atributo categórico todos los registros son "AllPub", excepto por uno "NoSeWa" y 2 NA. Dado que la casa con 'NoSewa' está en el juego de datos de entrenamiento, este atributo no ayudará en el modelado predictivo. Por tanto, podemos eliminarlo sin problema.
In [22]:
all_data = all_data.drop(['Utilities'], axis=1)
Functional: data description says NA means typical
In [23]:
all_data["Functional"] = all_data["Functional"].fillna("Typ")
Electrical : Tiene un valor NA. Pero este atributo en su mayoría es 'SBrkr', de manera que podemos reemplazar el valor perdido por 'SBrkr'.
In [24]:
all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])
KitchenQual: Sólo un valor NA, lo sustituimos por 'TA' (el más frecuente).
In [25]:
all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])
Exterior1st y Exterior2nd : De nuevo, Exterior 1 & 2 tienen un sólo valor perdido. Sustituimos por el más común.
In [26]:
all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])
SaleType : Rellenamos con el más frecuente ("WD")
In [27]:
all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])
MSSubClass : Na muy probablemente indica clase "no building", por lo que reemplazamos con None
In [28]:
all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")
Quedan aún valores perdidos?
In [29]:
# Verificamos si permanecen aún valores perdidos
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()
Out[29]:
Ya no quedan valores perdidos.
In [30]:
# MSSubClass = El tipo de edificios
all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)
# OverallCond
all_data['OverallCond'] = all_data['OverallCond'].astype(str)
# Año y mes de venta.
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)
In [31]:
from sklearn.preprocessing import LabelEncoder
cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond',
'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1',
'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',
'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond',
'YrSold', 'MoSold')
# procesa columnas, applicando LabelEncoder a los atributos categóricos
for c in cols:
lbl = LabelEncoder()
lbl.fit(list(all_data[c].values))
all_data[c] = lbl.transform(list(all_data[c].values))
# Shape
print('Shape all_data: {}'.format(all_data.shape))
In [32]:
# Adicionamos el total de pies cuadrados (TotalSF) de la vivienda
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']
In [33]:
numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index
# Verificamos el sesgo de todos los atributos numéricos
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\nSkew in numerical features: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(20)
Out[33]:
In [34]:
skewness = skewness[abs(skewness) > 0.75]
print("There are {} skewed numerical features to Box Cox transform".format(skewness.shape[0]))
from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
#all_data[feat] += 1
all_data[feat] = boxcox1p(all_data[feat], lam)
# all_data[skewed_features] = np.log1p(all_data[skewed_features])
In [35]:
all_data = pd.get_dummies(all_data)
print(all_data.shape)
In [36]:
train = all_data[:ntrain]
test = all_data[ntrain:]
In [37]:
# Mapa para detectar qué atributos están correlacionados con el precio de venta
corrmat = train.corr()
plt.subplots(figsize=(12,9))
sns.heatmap(corrmat, vmax=0.9, square=True)
Out[37]:
In [38]:
from sklearn.linear_model import LinearRegression, ElasticNet, Lasso, BayesianRidge, LassoLarsIC
from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
In [39]:
# Función de validación cruzada con mezcla aleatoria previa
n_folds = 5
def rmsle_cv(model):
kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)
rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))
return(rmse)
Regresión lineal (clásica):
In [40]:
REG = make_pipeline(RobustScaler(), LinearRegression())
Regresión LASSO:
In [41]:
LASSO = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))
Regresión Elastic Net:
In [42]:
ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))
Regresión con Kernel Ridge:
In [43]:
KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
Veamos ahora cómo estos modelos se desempeñan sobre los datos evaluando la métrica RMSLE.
In [44]:
score = rmsle_cv(REG)
print("Linear regression score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(LASSO)
print("Lasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(ENet)
print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(KRR)
print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
In [45]:
# Ajustamos el modelo con el mejor score a los datos y predecimos los valores
LASSO.fit(train, y_train)
y_pred = LASSO.predict(test)