Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Um exercício de Machine Learning usando a base "House Prices: Advanced Regression Techniques" do Kaggle

+1 voto
73 visitas
perguntada Dez 15, 2020 em Aprendizagem de Máquinas por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Este exercício é baseado na competição "House Prices: Advanced Regression Techniques" do Kaggle, em que o objetivo é prever o preço final (target) de uma casa na cidade de Ames, Iowa, com base nas vendas ocorridas entre Janeiro de 2006 e Julho de 2010, a partir de 79 variáveis explanatórias (features), tanto categóricas quando numéricas. É um problema de aprendizagem supervisionada, em que desejamos, baseado em um conjunto de dados de treinamento, que determinado modelo possa aprender a relação entre as features e o target, e assim ser capaz de prever o preço de venda de casas nesta cidade para novos dados. E, como a variável a ser prevista é contínua, nos deparamos com um problema de regressão.

Apesar de parte do nome da competição ser um tanto intimidador - Advanced Regression Techniques - na realidade a base disponível é um ótimo ponto de partida para quem está iniciando na área de aprendizado de máquina, pois apresenta desafios típicos em problemas de regressão, um rico conjunto de dados, e uma ótima documentação, disponíveis em www.kaggle.com/c/house-prices-advanced-regression-techniques

Compartilhe

6 Respostas

0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Análise Exploratória dos Dados (AED) - Parte 1

Nesta parte inicial exploraremos os dados disponíveis, buscando alguns insights que podem ser úteis para a etapa de modelagem. Desenvolvemos, nessa e nas demais etapas, notebooks no próprio kaggle, sob a competição House Prices: Advanced Regression Techniques, que serão detalhados ao longo dessa questão.

Setup e checagem inicial

Carregamos todas as bibliotecas que serão utilizadas ao longo dessa etapa, bem como setamos o padrão dos gráficos para o default do seaborn:

import os
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
from scipy.stats import norm

sns.set()
warnings.filterwarnings('ignore')

single_plot_w, single_plot_h = 1.5 * 3.8, 3.8
ylim = (0, 10 ** 6)

"Espiamos" o diretório da base de dados para verificar o que há disponível:

for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))

Obtendo:

/kaggle/input/house-prices-advanced-regression-techniques/sample_submission.csv
/kaggle/input/house-prices-advanced-regression-techniques/data_description.txt
/kaggle/input/house-prices-advanced-regression-techniques/train.csv
/kaggle/input/house-prices-advanced-regression-techniques/test.csv

Conforme descrição na própria página da base de dados, em train.csv temos os dados nominados, mas não há preços listados em test.csv, já que é o arquivo para submissão de previsões para quem participar da competição. Vamos então nos ater a train.csv, que será considerada como nossa única base neste exercício. Outro arquivo importante é data_description.txt, que contém a descrição de todos os features, como veremos a seguir. Carregamos então train.csv em um DataFrame, listamos suas colunas, e damos uma breve espiada em parte dos dados e em suas dimensões:

df = pd.read_csv("/kaggle/input/house-prices-advanced-regression-techniques/train.csv")

print("colunas:", df.columns)
print(df.head())

Que nos retorna:

colunas: Index(['Id', 'MSSubClass', 'MSZoning', 'LotFrontage', 'LotArea', 'Street',
       'Alley', 'LotShape', 'LandContour', 'Utilities', 'LotConfig',
       'LandSlope', 'Neighborhood', 'Condition1', 'Condition2', 'BldgType',
       'HouseStyle', 'OverallQual', 'OverallCond', 'YearBuilt', 'YearRemodAdd',
       'RoofStyle', 'RoofMatl', 'Exterior1st', 'Exterior2nd', 'MasVnrType',
       'MasVnrArea', 'ExterQual', 'ExterCond', 'Foundation', 'BsmtQual',
       'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinSF1',
       'BsmtFinType2', 'BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF', 'Heating',
       'HeatingQC', 'CentralAir', 'Electrical', '1stFlrSF', '2ndFlrSF',
       'LowQualFinSF', 'GrLivArea', 'BsmtFullBath', 'BsmtHalfBath', 'FullBath',
       'HalfBath', 'BedroomAbvGr', 'KitchenAbvGr', 'KitchenQual',
       'TotRmsAbvGrd', 'Functional', 'Fireplaces', 'FireplaceQu', 'GarageType',
       'GarageYrBlt', 'GarageFinish', 'GarageCars', 'GarageArea', 'GarageQual',
       'GarageCond', 'PavedDrive', 'WoodDeckSF', 'OpenPorchSF',
       'EnclosedPorch', '3SsnPorch', 'ScreenPorch', 'PoolArea', 'PoolQC',
       'Fence', 'MiscFeature', 'MiscVal', 'MoSold', 'YrSold', 'SaleType',
       'SaleCondition', 'SalePrice'],
      dtype='object')
   Id  MSSubClass MSZoning  LotFrontage  LotArea Street Alley LotShape  \
0   1          60       RL         65.0     8450   Pave   NaN      Reg   
1   2          20       RL         80.0     9600   Pave   NaN      Reg   
2   3          60       RL         68.0    11250   Pave   NaN      IR1   
3   4          70       RL         60.0     9550   Pave   NaN      IR1   
4   5          60       RL         84.0    14260   Pave   NaN      IR1   

  LandContour Utilities  ... PoolArea PoolQC Fence MiscFeature MiscVal MoSold  \
0         Lvl    AllPub  ...        0    NaN   NaN         NaN       0      2   
1         Lvl    AllPub  ...        0    NaN   NaN         NaN       0      5   
2         Lvl    AllPub  ...        0    NaN   NaN         NaN       0      9   
3         Lvl    AllPub  ...        0    NaN   NaN         NaN       0      2   
4         Lvl    AllPub  ...        0    NaN   NaN         NaN       0     12   

  YrSold  SaleType  SaleCondition  SalePrice  
0   2008        WD         Normal     208500  
1   2007        WD         Normal     181500  
2   2008        WD         Normal     223500  
3   2006        WD        Abnorml     140000  
4   2008        WD         Normal     250000  

[5 rows x 81 columns]

Podemos perceber que as colunas do dataframe seguem uma estrutura típica deste tipo de base de dados: a primeira coluna traz o 'Id' da observação, a última o nosso target, e entre elas 79 colunas de features, tanto numéricos quanto categóricos, que vamos explorar em seguida.

Já aproveitamos para verificar quais colunas possuem valores faltantes:

df.columns[df.isna().any()]

Obtendo:

Index(['LotFrontage', 'Alley', 'MasVnrType', 'MasVnrArea', 'BsmtQual',
       'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2',
       'Electrical', 'FireplaceQu', 'GarageType', 'GarageYrBlt',
       'GarageFinish', 'GarageQual', 'GarageCond', 'PoolQC', 'Fence',
       'MiscFeature'],
      dtype='object')

Ou seja, há uma série de features com valores faltantes, para os quais buscaremos um tratamento durante a etapa de modelagem. Partimos então para avaliação do target e das features.

Avaliação do target

Plotamos a distribuição do preço de venda das casas e verificamos medidas de assimetria e de curtose:

fig, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=(10, 6))
_ = sns.boxplot(df['SalePrice'], ax=ax_box)
_ = sns.distplot(df['SalePrice'], ax=ax_hist)
ax_box.set(xlabel='')
print("{} {:f}{} {} {:f}".format("Medida de assimetria (Skewness):", df['SalePrice'].skew(), ";", "Medida de curtose:", df['SalePrice'].kurt()))

Medida de assimetria (Skewness): 1.882876; Medida de curtose: 6.536282

A imagem será apresentada aqui.

Assim podemos verificar que nossa distribuição é Leptocúrtica Positivamente Assimétrica (medida de curtose acima de zero é relacionada com o afunilamento no pico em relação a Normal, e a medida de assimetria/skewness positiva indica o que verificamos visualmente: a cauda mais longa a direita). Interessante notar que tal distribuição provavelmente é a que esperaríamos em se tratando de preços de casas em uma cidade típica, em que podemos identificar claramente outliers / preços altos que se destoam dos demais.

Aplicamos então uma transformação logarítmica em nosso target, e verificamos que com isso basicamente a aproximamos para uma distribuição Normal, o que utilizaremos na etapa de modelagem:

fig, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=(10, 6))
logSalePrice = np.log(df['SalePrice'])
_ = sns.boxplot(logSalePrice, ax=ax_box)
_ = sns.distplot(logSalePrice, ax=ax_hist, fit=norm)
ax_box.set(xlabel='')
print("{} {:f}{} {} {:f}".format("Medida de assimetria (Skewness):", logSalePrice.skew(), ";", "Medida de curtose:", logSalePrice.kurt()))

Medida de assimetria (Skewness): 0.121335; Medida de curtose: 0.809532

A imagem será apresentada aqui.

Avaliação dos Features

Se verificarmos a base de dados apenas pelo tipo em que os dados são codificados:

df.dtypes.value_counts()

Obtemos:

object     43
int64      35
float64     3
dtype: int64

Mas não conseguimos na realidade concluir quais features são numéricas ou categóricas, pois a realidade da base é um pouco mais complexa, o que concluímos ao verificar a descrição das features disponível no arquivo data_description.txt.

Após analisar todas as features listadas, basicamente, podemos dividi-las em 18 features numéricos contínuos que se referem a dimensões / área (em pés / pés quadrados) (18 no total), 1 feature numéricos contínuos de valores monetários, 9 features numéricos discretos de quantidades referentes a alguma característica específica da casa, 5 features numéricos discretos relacionados a datas (com a observação que a feature 'GarageYrBlt' está codificada em float), 22 features categóricos ordinais (sendo que alguns estão codificados em valores inteiros e outros com a apenas a descrição) e 24 features categóricos nominais (sendo que um deles, MSSubClass, está codificado com valores inteiros, mas que não apresenta nenhuma ordinalidade).

Abaixo listamos os features em cada um desses grupos, com uma breve descrição:

Features numéricos contínuos que se referem a dimensões / área (em pés / pés quadrados) (18):
- 'LotFrontage': Pés lineares da rua conectados à propriedade;
- 'LotArea': Área do lote;
- 'MasVnrArea': Área de alvenaria;
- 'BsmtFinSF1': Área acabada do porão (do primeiro tipo);
- 'BsmtFinSF2': Área acabada do porão (do segundo tipo, se mais de um);
- 'BsmtUnfSF': Área não acabada do porão;
- 'TotalBsmtSF': Área total do porão;
- '1stFlrSF': Área do primeiro andar;
- '2ndFlrSF': Área do segundo andar;
- 'LowQualFinSF': Área com acabamento de baixa qualidade (todos os andares);
- 'GrLivArea': Área útil acima do nível solo;
- 'GarageArea': Área da garagem;
- 'WoodDeckSF': Área do deck de madeira;
- 'OpenPorchSF': Área de varanda aberta;
- 'EnclosedPorch': Área de varanda fechada;
- '3SsnPorch': Área da varanda de três temporadas;
- 'ScreenPorch': Área da tela da varanda;
- 'PoolArea': Área da piscina;

Features numéricos contínuos de valores monetários (1):
- 'MiscVal': Valor monetário de características diversas não incluídas em outras categorias;

Features numéricos discretos de quantidades de alguma característica da casa (9):
- 'BsmtFullBath': Banheiros completos no porão;
- 'BsmtHalfBath': Lavabos no porão;
- 'FullBath': Banheiros completos acima do nível do solo;
- 'HalfBath': Lavabos completos acima do nível do solo;
- 'BedroomAbvGr': Quartos acima do nível do solo (não inclui quartos no subsolo);
- 'KitchenAbvGr': Cozinhas acima do nível do solo;
- 'TotRmsAbvGrd': Total de cômodos acima do nível do solo (não inclui banheiros);
- 'Fireplaces': Número de lareiras;
- 'GarageCars': Capacidade da garagem (número de carros);

Features numéricos discretos relacionados a datas (5) (exceção para 'GarageYrBlt', que possui valores contínuos):
- 'YearBuilt': Ano de construção da casa;
- 'YearRemodAdd': Ano da reforma da casa (mesma do ano da construção de não houver);
- 'GarageYrBlt': Ano de construção da garagem;
- 'MoSold': Mês de venda;
- 'YrSold': Ano de venda;

Features categóricos ordinais (22):
- 'LotShape': Formato geral do lote (regular, ligeiramente irregular, moderadamente irregular, irregular);
- 'Utilities': Utilidades disponíveis (de apenas eletricidade a todas);
- 'LandSlope': Inclinação do terreno (3 níveis);
- 'OverallQual': Qualidade geral da casa (de 0 a 10);
- 'OverallCond': Condições gerais da cada (de 0 a 10);
- 'ExterQual': Qualidade do material no exterior (5 níveis);
- 'ExterCond': Avaliação a condição atual do material no exterior (5 níveis);
- 'BsmtQual': Avaliação a altura do porão (5 níveis);
- 'BsmtCond': Avaliação do estado geral do porão (5 níveis);
- 'BsmtFinType1': Classificação da área acabada do porão (6 níveis + NA);
- 'BsmtFinType2': Classificação da área acabada do porão (para o segundo tipo, se mais de um) (6 níveis + NA);
- 'BsmtExposure': Avaliação da exposição externa do porão;
- 'HeatingQC': Qualidade e condição do aquecimento (5 níveis);
- 'KitchenQual': Qualidade da cozinha (5 níveis);
- 'Functional': Funcionalidade inicial (presuma típica, a menos que as deduções sejam garantidas) (8 níveis);
- 'FireplaceQu': Qualidade da lareira (5 níveis + NA);
- 'GarageFinish': Estado do acabamento interior da garagem (acabado, mal acabado, não acabado, sem garagem, NA);
- 'GarageQual': Qualidade da garagem (5 níveis + NA);
- 'GarageCond': Condição da garagem (5 níveis + NA);
- 'PavedDrive': Pavimento da pista até a garagem (pavimentado, parcialmente pavimentado, sem pavimento/cascalho);
- 'PoolQC': Qualidade da piscina (4 níveis + NA);
- 'Fence': Qualidade da cerca (4 níveis + NA);

Features categóricos nominais (24):
- 'Street': Tipo de estrada de acesso à propriedade (pavimentada ou de cascalho);
- 'Alley': Tipo de acesso à propriedade (pavimentada, de cascalho ou NA);
- 'MSSubClass': Classificação do tipo de imóvel;
- 'MSZoning': Classificação geral do zoneamento (Comercial, Agricultural, Residencial de média densidade, etc);
- 'LandContour': Tipo de nivelamento do terreno da propriedade;
- 'LotConfig': Configuração do lote (ex: esquina, dentro do lote, etc.);
- 'Neighborhood': Bairro;
- 'Condition1': Proximidade a alguma via específica (ex: Adjacente à Ferrovia Norte-Sul);
- 'Condition2': Idem a anterior, se mais de uma;
- 'BldgType': Tipo de moradia (ex: duplex, unidade interna da casa geminada, etc.);
- 'HouseStyle': Estilo de moradia (ex: casa térrea, sobrado, etc.). Semelhante a MSSubClass, mas menos detalhado;
- 'RoofStyle': Tipo de telhado;
- 'RoofMatl': Material do telhado;
- 'Exterior1st': Tipo de acabamento exterior;
- 'Exterior2nd': Idem a anterior, se mais de um;
- 'MasVnrType': Tipo da cobertura da alvenaria;
- 'Foundation': Tipo de fundação;
- 'Heating': Tipo de aquecimento;
- 'CentralAir': Ar condicionado central (sim ou não);
- 'Electrical': Sistema elétrico;
- 'GarageType': Localização da garagem em relação a casa;
- 'MiscFeature': Características diversas não incluídas em outras categorias;
- 'SaleType': Tipo de venda (tipo de garantia, tipo de contrato);
- 'SaleCondition': Condição de venda (ex: normal, venda entre membros de família, etc.).

Análise dos Features e de sua relação com o target

Nessa etapa plotamos as distribuições de cada features, bem como os gráficos de dispersão, em relação ao preço de venda, para as features numéricas, e bloxplots para as features categóricas, divididos conforme grupos citados anteriormente.

Assim, para os os features numéricos que se referem a dimensões / área:

def plot_numeric(feature_list, ylim=(0, 10 ** 6), single_plot_w=1.5 * 3.8, single_plot_h=3.8):
    m, n = len(feature_list), 2
    fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w, m * single_plot_h))
    for i, feature in enumerate(feature_list):
        g = sns.distplot(df[feature], ax=ax[i][0], kde=False)
        g = sns.regplot(x=df[feature], y=df['SalePrice'], ax=ax[i][1])
        g.set(ylim=ylim)

area_features = ['LotFrontage', 'LotArea', 'MasVnrArea', 'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF', '1stFlrSF', '2ndFlrSF', 'LowQualFinSF', 'GrLivArea', 'GarageArea', 'WoodDeckSF', 'OpenPorchSF', 'EnclosedPorch', '3SsnPorch', 'ScreenPorch', 'PoolArea']        
plot_numeric(area_features)

O que nos gera o seguinte gráfico:

Features numéricos que se referem a dimensões / área

Separamos aqui dois exemplos interessantes:

A imagem será apresentada aqui.

Podemos traçar algumas conclusões com relação a estes gráficos:

  • Como esperado, o preço de venda possui relação aproximadamente linear com a área útil total acima do nível do solo (GrLivArea), mas com considerável heteroscedasticidade. Também podemos observar alguns outliers (valores para áreas acima de 4000 pés quadrados), mas que segundo [ 1 ] não representam observações incorretas, e sim valores não usuais de venda;

  • Outros features de área também possuem relação aproximadamente linear com o preço de venda, por exemplo, área total do porão;

  • Moradores de Ames parecem não ser muito fãs de piscinas...

Observação: notamos que, se utilizarmos a transformação logarítmica em SalePrices, reduzimos também heteroscedasticidade em relação às medidas de área. Por exemplo:

fig, ax = plt.subplots(figsize=(single_plot_w, single_plot_h))
g = sns.regplot(x=df['GrLivArea'], y=np.log(df['SalePrice']))

A imagem será apresentada aqui.

Com relação a valores monetários, possuímos apenas um feature, correspondente ao valor de itens adicionais que uma casa possuí, o qual plotamos por:

m, n = 1, 2
fig, (ax1, ax2) = plt.subplots(m, n, figsize=(n * single_plot_w, m * single_plot_h))
g1 = sns.distplot(df['MiscVal'], ax=ax1, kde=False)
g2 = sns.regplot(x=df['MiscVal'], y=df['SalePrice'], ax=ax2)
_ = g2.set(ylim=ylim)

A imagem será apresentada aqui.

Aparentemente, os valores de itens adicionais, individualmente, parecem pouco informativos sobre o valor do imóvel.

Para os features numéricos discretos de quantidades de itens específicos da casa:

def plot_categorical(feature_list, single_plot_w=1.5 * 3.8, single_plot_h=3.8):
    m, n = len(feature_list), 2
    fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w, m * single_plot_h))
    for i, feature in enumerate(feature_list):
        g = sns.countplot(df[feature], ax=ax[i][0])
        g = sns.boxplot(x=df[feature], y=df['SalePrice'], ax=ax[i][1])       
        g.set_yticklabels(['{:,.0f}'.format(y) + 'K' for y in g.get_yticks()/1000])

discrete_features = ['BsmtFullBath', 'BsmtHalfBath', 'FullBath', 'HalfBath', 'BedroomAbvGr', 'KitchenAbvGr', 'TotRmsAbvGrd', 'Fireplaces', 'GarageCars']
plot_categorical(discrete_features)

O que nos gera o seguinte plot:

Features numéricos discretos de quantidades de ítens específicos da casa

Em que separamos dois exemplos interessantes:

A imagem será apresentada aqui.

Podemos verificar que algumas quantidades possuem boa correlação com o preço do imóvel (ex: número de banheiros, quando diferente de zero), apesar de, em alguns casos, a relação para determinados valores não ser óbvia, possivelmente pela interação com outra variável ou pelo tamanho da amostra. Por exemplo, o número de quartos acima do nível do solo, em média, é correlacionado com o preço do imóvel quando de dois a quatro quartos, mas tal correlação não se mantém para o próximo quarto. Porém, o número de amostras para 1 ou 5 quartos é muito menor em relação às de 2 a 4 quartos. Algo semelhante ocorre para cômodos totais acima do nível do solo: temos uma relação quase linear em média, até chegarmos a 12 ou mais cômodos, quando também o número de amostras é bem menor. Para garagens, o mesmo, até a terceira garagem, mas uma quarta "reduz" o preço do imóvel (mas neste caso quase não há observações). Outro fato interessante é que ter duas cozinhas não parece ser um bom negócio...

0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Análise Exploratória dos Dados (AED) - Parte 2

Para os features relacionados a datas:

date_features = ['YearBuilt', 'YearRemodAdd', 'GarageYrBlt', 'MoSold', 'YrSold']

m, n = 3, 1

fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w * 4, m * single_plot_h * 2))
for feature, subplot in zip(date_features[:3], ax.flatten()):
    g = sns.countplot(df[feature], ax=subplot)
    g.set_xticklabels(g.get_xticklabels(), rotation=90)

Features relacionados a datas (distribuição)

fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w * 4, m * single_plot_h * 2))
for feature, subplot in zip(date_features[:3], ax.flatten()):
    g = sns.boxplot(x=df[feature], y=df['SalePrice'], ax=subplot)
    g.set_xticklabels(g.get_xticklabels(), rotation=90)
    g.set_yticklabels(['{:,.0f}'.format(y) + 'K' for y in g.get_yticks()/1000])

Features relacionados a datas (boxplot)

plot_categorical(date_features[3:])

A imagem será apresentada aqui.

Como possíveis conclusões:

  • Todas as casas construídas antes de 1950, que não tiveram reforma após isso, estão com ano da última reforma em 1950 (lembrando que, se não houve reforma, esse ano deveria indicar o ano de construção, conforme descrição da documentação).

  • É possível notar um vácuo de construções na década de 80...

  • Construções recentes tendem a sinalizar, em média, um maior preço de venda, conforme seria realmente esperado.

  • O ano de venda não parece ser muito informativo sobre o preço de venda.

  • O mercado é mais aquecido no verão americano.

Já para os features ordinais categóricos, antes de plotar, estabelecemos a ordenação adequada para cada um (infelizmente, cada feature, em geral, possui uma ordenação diferente do outro), com exceção para alguns de qualidade:

category_ordinal_features = ['LotShape', 'Utilities', 'LandSlope', 'OverallQual', 'OverallCond', 'ExterQual', 'ExterCond', 'BsmtQual', 'BsmtCond', 'BsmtFinType1', 'BsmtFinType2', 'BsmtExposure', 'HeatingQC', 'KitchenQual', 'FireplaceQu', 'GarageFinish', 'GarageQual', 'GarageCond', 'PoolQC', 'Fence']
def plot_categorial_ordered(feature_list, order_list, single_plot_w=1.5 * 3.8, single_plot_h=3.8):
    m, n = len(feature_list), 2
    fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w, m * single_plot_h))
    for i, feature in enumerate(feature_list):
        if order_list[i] != "":
            g = sns.countplot(df[feature], ax=ax[i][0], order=order_list[i])
            g = sns.boxplot(x=df[feature], y=df['SalePrice'], ax=ax[i][1], order=order_list[i])
            _ = g.set_yticklabels(['{:,.0f}'.format(y) + 'K' for y in g.get_yticks()/1000])
        else:
            g = sns.countplot(df[feature], ax=ax[i][0])
            g = sns.boxplot(x=df[feature], y=df['SalePrice'], ax=ax[i][1])
            _ = g.set_yticklabels(['{:,.0f}'.format(y) + 'K' for y in g.get_yticks()/1000])

order_list = []
for feature in category_ordinal_features:
    if feature in ['ExterQual', 'PoolQC']:
        order_list.append(["Fa", "TA", "Gd", "Ex"])
    elif feature in ['BsmtQual', 'BsmtCond', 'HeatingQC', 'KitchenQual', 'GarageQual', 'GarageCond', 'FireplaceQu', 'ExterCond']:
        order_list.append(["Po", "Fa", "TA", "Gd", "Ex"])
    elif feature in ['BsmtExposure']:
        order_list.append(["No", "Mn", "Av", "Gd"])
    elif feature in ['BsmtFinType1', 'BsmtFinType2']:
        order_list.append(["Unf", "LwQ", "Rec", "BLQ", "ALQ", "GLQ"])
    elif feature in ['GarageFinish']:
        order_list.append(["Unf", "RFn", "Fin"])
    elif feature in ['PavedDrive']:
        order_list.append(["N", "P", "Y"])
    elif feature in ['Fence']:
        order_list.append(["MnWw", "GdWo", "MnPrv", "GdPrv"])
    else:
        order_list.append("")

plot_categorial_ordered(category_ordinal_features, order_list)

Que nos gera:

Features categóricos ordinais

E como um subconjunto interessante:

A imagem será apresentada aqui.

Algumas conclusões:

  • Alguns features de qualidade possuem interessante relação (aproximadamente linear ou exponencial) com o preço médio, como condições do porão, acabamento da garagem, qualidade geral da casa, qualidade do porão, condições do porão, qualidade do exterior e qualidade da cozinha, mas em geral com visível heteroscedasticidade.

  • Uma surpresa é o feature de condições gerais da casa, que apresenta os maiores preços médios de venda para valores iguais a 5 ou 9, (e sem diferenças significativas entre 6 e 8). Notar que a maioria das observações se concentra no valor 5.

  • Algumas medidas de qualidade / condições são poucas informativas diretamente sobre preços para faixas de classificação de mais de um valor (ex: cercas: só há diferenciação para a última medida de qualidade, porém com poucas observações).

  • Alguns features possuem quase a totalidade das observações concentradas em apenas um valor: quase todas as casas tem todas as utilidades previstas na pesquisa, ou quase todas as garagens são avaliadas como TA (Average/Typical), tanto em relação a qualidade quanto as condições atuais. O mesmo ocorre na avaliação das condições do porão e nas condições externas da cada. Ou seja, as classificação de qualidade para "condições", que apresentam o campo TA (Average/Typical) como possível resposta, tendem a concentrar todas as notas neste valor, talvez por ser conveniente, em termos de avaliação de "condições", responder que o mesmo está em condições médias / típicas na avaliação do imóvel. Mas isso não ocorre na avaliação de qualidade do material (original) de qual o item é construído.

Para os features categóricos nominais, teremos:

category_nominal_features = ['Neighborhood', 'Exterior1st', 'Exterior2nd', 'MSSubClass', 'MSZoning', 'Street', 'Alley', 'PavedDrive', 'LandContour', 'LotConfig', 'Condition1', 'Condition2', 'BldgType', 'HouseStyle', 'RoofStyle', 'RoofMatl', 'MasVnrType', 'Foundation', 'Heating', 'CentralAir', 'Electrical', 'GarageType', 'MiscFeature', 'Functional', 'SaleType', 'SaleCondition']

m, n = 3, 1

fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w * 2, m * single_plot_h * 1.5))
for feature, subplot in zip(category_nominal_features[:3], ax.flatten()):
    g = sns.countplot(df[feature], ax=subplot)
    g.set_xticklabels(g.get_xticklabels(), rotation=22.5)

Features categóricos nominais - grupo 1 - distribuições

fig, ax = plt.subplots(m, n, figsize=(n * single_plot_w * 2, m * single_plot_h * 1.5))
for feature, subplot in zip(category_nominal_features[:3], ax.flatten()):
    g = sns.boxplot(x=df[feature], y=df['SalePrice'], ax=subplot)
    g.set_yticklabels(['{:,.0f}'.format(y) + 'K' for y in g.get_yticks()/1000])
    g.set_xticklabels(g.get_xticklabels(), rotation=22.5)

Features categóricos nominais - grupo 1 - boxplot

plot_categorical(category_nominal_features[3:])

Features categóricos nominais - grupo 2 - distro e boxplot

Separando o exemplo dos bairros:

A imagem será apresentada aqui.
A imagem será apresentada aqui.

Algumas características claramente são capazes de diferenciar diretamente o preço médio das casas. Por exemplo:

  • Bairros parecem ser interessante preditores de preços de casas, mas alguns bairros possuem poucas amostras.

  • Material de cobertura exterior da casa, se dos tipos Metal Siding e Wood Siding, que representam razoável parte da amostra, tendem a caracterizar casas mais baratas em relação as que utilizam Vinyl Siding (que é a maioria).

  • O tipo da área em que se encontram são bons preditores de preço médio: casas em áreas comerciais possuem preços menores, enquanto casas em vilarejos especiais (Floating Village) em geral são mais caras as demais.

  • A maioria das casas possui um entre dois tipos de fundação (Brick & Tile ou Cinder Block), sendo que a variância de preços relativa a cada um é baixa e os preços médios são bem distintos entre elas.

  • Não possuir ar central é um bom indicativo de que a casa será pouco valorizada.

  • Tipo da cobertura de alvenaria, quando de pedra, indica valorização no preço na casa;

  • Casas com garagens não anexadas na casa tendem a ser mais baratas;

  • Casas com sistema elétrico atualizado (padrão) tendem a valer mais;

  • Casas vendidas na modalidade Partial (não estavam concluídas no momento da venda, logo podemos associar com casas novas) são vendidas a um preço em geral superior as demais, como esperado.

  • A maioria das casas é da categoria 20 ('MSSubClass'), que correspondem a casas de um andar apenas, construídas a partir de 1946. E MSSubClass parece ser um bom preditor para os preços.

  • Casas em que o acesso é por estrada pavimentada possuem maior preço médio.

Em alguns features categóricos, determinada categoria acaba sendo extremamente dominante, como no caso de uma segunda característica de proximidade a alguma via importante (Condition2), material do telhado ou sistema de aquecimento.

0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Análise Exploratória dos Dados (AED) - Parte 3

Análise das correlações entre os features numéricos

Uma forma interessante de verificar as correlações entre as features e entre elas e o target é por meio do mapa de calor. Assim, plotamos o mapa de calor completo das variáveis numéricas (e categóricas ordinais já representadas numericamente) por:

corr_matrix = df.drop(['Id', 'MSSubClass'], axis=1).corr()
fig, ax = plt.subplots(figsize=(24, 18))
_ = sns.heatmap(corr_matrix, annot=True, fmt='.1f', vmin=-1, vmax=1, center= 0)

Que nos gera o mega mapa de calor:

Mapa de calor das correlações entre os features numéricos

Apenas para exemplificar, em uma versão reduzida do mapa de calor, plotamos, por exemplo, apenas as features que possuem correlação maior que 0.5 com SalePrice:

fig, ax = plt.subplots(figsize=(9, 6))
threshold = .5
sns.set(font_scale=.8)
corr_matrix = df.corr()
_ = sns.heatmap(df[corr_matrix[corr_matrix['SalePrice'] > threshold]['SalePrice'].index.tolist()].corr(), annot=True, fmt='.2f', vmin=-1, vmax=1, center= 0, cbar=False)

A imagem será apresentada aqui.

Algumas conclusões:

  • Os features com as maiores correlações com o preço de venda são qualidade geral e área útil acima do solo, mas também merecem atenção área do porão, número de banheiros, área / capacidade da garagem, ano de construção, tipo de cobertura da alvenaria e quantidade de quartos, como é razoável de se supor se imaginamos o que procuramos em uma casa. Um fato curioso é a correlação entre o número de lareiras e o preço do imóvel. Já a correlação entre a área do terreno e a área útil acima do solo é apenas na ordem de 0.3.

  • Assim como não são correlacionados com o preço, como vimos anteriormente, algumas features quantitativas relacionadas a varandas, piscina, valor de itens diversos, mês e ano da venda também possuem baixíssima ou nenhuma correlação com os demais features ('3SsnPorch', 'ScreenPorch', 'PoolArea', 'MiscVal', 'MsSold', 'YrSold'). O mesmo vale para 'BsmtFinSF2' (segunda área acabada do porão), 'LowQualFinSF' (área útil de baixa qualidade), 'BsmtHalfBath' (Lavabos)

  • Notamos alta correlação entre alguns pares de features, em que o par possui capacidade similar de previsão de valor de venda, indicativo de multicolinearidade. São exemplos a área da garagem e sua capacidade, ano da construção da casa e ano de construção da garagem, área total útil e número de cômodos acima do solo, área do primeiro andar e área do porão (visto que para estes pares de features, em geral, uma é consequência natural da outra).

  • A condição geral da casa não é correlacionada com o preço devido ao fato que já apontamos (maioria das observações concentrada na nota 5, e que apresentam maior preço médio que as notas vizinhas, gerando um padrão não linear). Porém a condição geral é negativamente correlacionada com o ano da construção, o que seria um fato interessante a se supor.

  • Apenas por curiosidade, quanto mais nova a casa, menor são as áreas destinadas a varandas fechadas...

Porém, o mapa de calor só é capaz de indicar multicolinearidades entre duas features, e não entre combinações lineares. Assim, como identificar multicolinaridade perfeita não presente no mapa de calor, se existir? Com base nas features disponíveis, é de se supor que:

BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF = TotalBsmtSF (área total do porão é igual as áreas de partes específicas deste presentes na base)

1stFlrSF + 2ndFlrSF + LowQualFinSF = GrLivArea (área total acima do solo é igual a soma da área dos andares e da área declarada como não acabada, se supormos que, para casas com mais de dois andares, as áreas dos andares superiores são imputadas na área do segundo andar...)

Podemos mostrar que sim, isto ocorre, por:

print(df['TotalBsmtSF'].corr(df['BsmtFinSF1'] + df['BsmtFinSF2'] + df['BsmtUnfSF']))
print(df['GrLivArea'].corr(df['1stFlrSF'] + df['2ndFlrSF'] + df['LowQualFinSF']))

1.0
1.0

Um gráfico interessante é o de dispersão entre as features, conforme apresentado abaixo (excluindo-se outliers em que área útil acima do solo é maior que 4000 pés quadrados, apenas para melhorar a visualização das relações):

drop_list=['Id', 'MSSubClass', 'BsmtFinSF2', 'BsmtUnfSF', '2ndFlrSF', 'LowQualFinSF', '3SsnPorch', 'ScreenPorch', 'PoolArea', 'MiscVal', 'MoSold', 'YrSold', 'BsmtFinSF2', 'LowQualFinSF', 'BsmtHalfBath', 'KitchenAbvGr']
g = sns.pairplot(df[df['GrLivArea'] < 4000].drop(drop_list, axis=1))

Mega plot de dispersão entre features numéricas

PS: a melhor forma de visualizar o gráfico acima é fazer o download da imagem original via Postimages.

E para ilustrar, o mesmo plot para algumas das features numéricas com maior correlação com o preço de venda:

n = 6
cols = corr_matrix.drop(['GarageCars']).nlargest(n, 'SalePrice').index.tolist()
if 'GrLivArea' in cols:
    sns.pairplot(df[df['GrLivArea'] < 4000][cols]) #, size=2.5);

A imagem será apresentada aqui.

Com relação a esses gráficos, alguns fatos interessantes:

  • Casas classificadas como de melhor qualidade tendem a ser as de maior área útil acima do solo e com porões e garagens maiores;
  • Fica visível como a área do porão é geralmente limitada pela área do primeiro andar (o que é esperado, claro), e há alguns raros curiosos outliers em que a área do porão é maior que a área útil do primeiro piso.

Análise das associações entre os features categóricos

Para a análise de associações entre os features categóricos e entre esses e o preço de venda utilizamos a biblioteca dython de Shaked Zychlinski. Basicamente a função associations() gera o mapa de calor em função do coeficiente de incerteza de Theil's U, que é assimétrico (x ~ y \(\neq\) y ~ x), o que permite explorar relações interessantes entre as features nominais, bem como a razão de correlação entre as features categóricas e 'SalePrices' (mais detalhes sobre esse tipo de abordagem pode ser verificado em um breve post de Zychlinski no towardsdatascience.com).

Assim, instalamos a biblioteca dython por (já que não está instalada por padrão no kaggle):

!pip install git+https://github.com/shakedzy/dython.git

E para a geração do mapa de calor para variáveis categóricas ordinais que não estão codificadas numericamente:

from dython import nominal

nominal.associations(df[category_ordinal_features + ['SalePrice']].drop(["OverallQual", "OverallCond"], axis=1), fmt='.1f', theil_u=True, figsize=(12, 9));

A imagem será apresentada aqui.

Podemos tirar algumas conclusões:

  • Qualidade Externa, do Porão e da Cozinha são os itens mais associados com o valor do imóvel;
  • O acabamento da garagem diz muito sobre sua qualidade e condições, mas o contrário não é verdade;
  • Apenas como curiosidade, quem se preocupa com a qualidade da cozinha também se preocupa com a qualidade externa, mas não muito com a do porão... E a qualidade da cerca nos diz alguma coisa sobre a qualidade da piscina, mas o inverso não é verdade, bem como o acabamento da garagem diz muito sobre sua qualidade e condições, mas o inverso não é verdade...

E para a geração do mapa de calor para variáveis categóricas nominais:

_ = nominal.associations(df[category_nominal_features + ['SalePrice']], theil_u=True, figsize=(15, 12), fmt='.1f', cbar=False)

Que nos gera:

A imagem será apresentada aqui.

Podemos verificar que alguns features possuem associações significativas em ambas as direções, por exemplo:

  • Bairro e a Classificação do tipo da casa;
  • Tipo e Estilo da Casa com Classificação do tipo da casa, sendo que o Tipo e Classificação podem ser até ser redundantes de certa forma;
  • Condição de Venda e Tipo de Venda;

Alguns features possuem associação mais significativa em uma direção, mas não na outra:

  • O bairro pode explicar parte do tipo de classificação da zona em que se encontra a casa (comercial, zona agrícola, residencial), mas pelo tipo de zoneamento não é possível identificar o bairro;
  • O bairro permite prever parte da fundação que será utilizada, mas o inverso não é verdade;
  • O estilo do telhado ('RoofStyle') permite prever o seu material ('RoofMatl'), mas o inverso não é verdade;

Alguns features aparentam ser bons preditores para o preço de venda (Tipo de Fundação, Tipo da Garagem, Tipo Cobertura da Alvenaria ('MasVnrType'), Cobertura Exterior da Casa ('Exterior1st e 2nd')).

Itens diversos ('MscFeature'), 2o tipo de condições de proximidade ('Condition2') e tipo da rua ('Street') não têm poder significativo de explicar alguma feature ou o target.

Podemos também gerar um mega mapa com todas as associações entre features categóricas ordinais e nominais, incluindo "OverallQual", "OverallCond", que são representadas numericamente, por:

_ = nominal.associations(df[category_ordinal_features + category_nominal_features], theil_u=True, fmt='.1f', figsize=(30 * .8, 24 * .8), cbar=False)

Que nos gera:

Mega mapa com todas as associações entre features categóricas ordinais e nominais

Conclusões:

  • A qualidade geral da casa se associa com várias features categóricas nominais, como Tipo da Cobertura da Alvenaria, Tipo da Fundação e Tipo da Garagem, mas principalmente com o bairro em que se localiza (pelo bairro, você pode esperar qual a qualidade da casa que vai encontrar);

  • O bairro explica muito das características da casa;

  • Além dessa, não há muita associação entre as features categóricas ordinais e nominais, com exceção do tipo de garagem e sua avaliação de qualidade / condições / acabamento.

Referências:

[ 1 ] Cock, Dean. Ames, Iowa: Alternative to the Boston Housing Data as an End of Semester Regression Project. Journal of Statistics Education, Volume 19, Number 3, 2011.

0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Engenharia de Features e Modelagem (Parte 1)

Iniciamos essa etapa carregando as bibliotecas que serão utilizadas e as configurações iniciais, além da carga da base de dados (train.csv).

Observação: trabalharemos aqui apenas com a base nomeada, realizando nossas avaliações por validação cruzada. Para submissão à competição, bastaria realizar as transformações aqui apresentadas ao conjunto train.csv + test.csv, treinar o modelo desejado sobre train, realizar as previsões sobre os dados de teste e submeter.

import os
import datetime
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.preprocessing import OrdinalEncoder
from sklearn.model_selection import KFold, train_test_split, GridSearchCV, RandomizedSearchCV, cross_validate
from sklearn.inspection import permutation_importance

from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.ensemble import GradientBoostingRegressor

# Dribla warnings do sklearn:
def warn(*args, **kwargs):
    pass
import warnings
warnings.warn = warn

sns.set()  # seta gráficos nas configurações padrões do seaborn
random_state = 0  # para resultados serem reproduzíveis
number_splits = 5  # número de divisões para o Kfold

df = pd.read_csv("/kaggle/input/house-prices-advanced-regression-techniques/train.csv")

As métricas utilizadas para avaliação dos modelos será a raiz quadrática média dos erros (RMSE) entre os valores observados e as predições bem como o \(R^2\) . Para fins de aprendizado, e também curiosidade, avaliaremos mais de um modelo industrial de regressão do sklearn, e para tanto criamos algumas funções para modelagem, avaliação e apresentação dos resultados de forma que a inclusão e exclusão de algum modelo possua um baixo custo:

score_list = ['neg_mean_squared_error','r2']
limits = [10.5, 14]


def model_string_formatter(string):
    return string.split("(", 1)[0]


def rmse(neg_mean_absolute_error):
    return np.sqrt(-neg_mean_absolute_error)


def plot_comparison(results):    

    g = sns.PairGrid(results, y_vars="Model", x_vars=["RMSE", "r2"], height=(len(results) * 0.8), aspect= (4 / len(results)))
    g.map(sns.barplot, palette="Blues_r")
    sns.set()
    for graph, ax in enumerate(g.axes.flat):
        for tick in range(len(results)):
            ax.text(results.iloc[:, graph + 1].values[tick], tick, '{:.3f}'.format(results.iloc[:, graph + 1].values[tick]), color='black', verticalalignment="center")
    plt.show()


def evaluate(models, X, y, kfold, scoring=score_list):

    rmse_perf, r2_perf, model_desc, cv_results = [], [], [], []

    for i, model in enumerate(models):

        cv_results.append(cross_validate(model, X, y, cv=kfold, scoring=score_list, return_train_score=True))

        rmse_perf.append(rmse(np.mean(cv_results[i]['test_neg_mean_squared_error'])))
        r2_perf.append(np.mean(cv_results[i]['test_r2']))
        model_desc.append(model_string_formatter(str(model)))

    results = pd.DataFrame({"Model": model_desc, "RMSE": rmse_perf, "r2": r2_perf, "model": models})
    results.sort_values(by=['RMSE'], ascending=True, inplace=True)
    results = results.reset_index(drop=True)
    plot_comparison(results)
    return results


def plot_prediction(X, y, results, kfold, model):

    chosen_one = results[results["Model"] == model]
    train_index, test_index = next(kfold.split(X, y))
    X_test, y_test = X.iloc[test_index], y.iloc[test_index]
    X_train, y_train = X.iloc[train_index], y.iloc[train_index]
    fig, ax = plt.subplots(figsize=(6,6))
    predictions = chosen_one["model"].values[0].fit(X_train, y_train).predict(X_test)
    _ = sns.regplot(y_test, y=predictions, ax=ax)
    ax.set(ylim=limits, xlim=limits, xlabel='Actual Price', ylabel='Predicted Price', title=model)
    plt.plot(limits, limits, color='r', linestyle='-', linewidth=1, alpha=.7)
    overlay = 'RMSE(mean): {:.3f}\nr2(mean): {:.3f}'.format(chosen_one['RMSE'].to_numpy()[0], chosen_one['r2'].to_numpy()[0])
    plt.annotate(s=overlay,xy=(12, 10.75) ,size='large')
    plt.show()

Estabelecendo uma referência

Definimos a seguinte lista de modelos, inicialmente sem customização dos hiperparâmetros (ou seja, vamos testá-los inicialmente com os hiperâmetros padrões, para cada um, do sklearn):

models = [
    LinearRegression(),
    Ridge(random_state=random_state),
    Lasso(random_state=random_state),
    GradientBoostingRegressor(random_state=random_state),
    ]

Criamos uma versão Mickey Mouse da base de dados, com apenas os features numéricos, excluíndo todos os features com algum dado faltante, apenas para realizarmos um teste e verificarmos qual seria o desempenho se nada mais fosse feito com a base de dados:

X = df.drop(['Id', 'SalePrice'], axis=1)
y = np.log(df['SalePrice'])

X.dropna(axis=1, inplace=True)
X = X.select_dtypes(exclude=['object'])

E assim testamos os modelos, com validação cruzada, em suas versões padrão, por meio das funções desenvolvidas acima, e obtemos as médias dos resultados, sobre os dados de teste, de cada divisão da validação cruzada (porém, o plot da previsão x valores atuais de preços é apenas para ilustração, e se refere a divisão dos dados em treinamento/testes do primeiro fold da validação cruzada):

kfold = KFold(n_splits=number_splits, shuffle=True, random_state=random_state)

results = evaluate(models, X, y, kfold)

plot_prediction(X, y, results, kfold, model="LinearRegression")

A imagem será apresentada aqui.
A imagem será apresentada aqui.

Podemos observar que o Gradient Boosting Regressor apresenta um resultado razoável sobre a versão Mickey Mouse de nossa base, sem qualquer parametrização customizada, não à toa parece ser o modelo industrial querido das competições de ML. O modelo de regressão linear está um pouco aquém, e o interessante é que é o Lasso apresenta um performance sofrível, o que podemos supor inicialmente que seja devido ao não tuning dos hiperparâmetros, o que vamos avaliar mais em frente (não vamos desistir dele).

Com relação ao campeão atual, fazemos um exercício de permutation importance, ou seja, qual seria a queda na métrica \(R^2\) pela permutação aleatória de determinado feature (quebrando-se assim a correlação real entre ele o target):

def importance_eval(model, X, y, kfold):

    train_index, test_index = next(kfold.split(X, y))
    X_train, y_train = X.iloc[train_index], y.iloc[train_index]
    X_test, y_test = X.iloc[test_index], y.iloc[test_index]
    n_repeats = 30

    importances = permutation_importance(model.fit(X_train, y_train), X_test, y_test, n_repeats=n_repeats, random_state=random_state)
    for i in importances.importances_mean.argsort()[::-1]:
        if (importances.importances_mean[i] - 2 * importances.importances_std[i] > 0) & (importances.importances_mean[i] > 0.005):
            print("{}: {:.3f} +/- {:.3f}".format(X_train.columns.values[i], importances.importances_mean[i], importances.importances_std[i]))

model = GradientBoostingRegressor(random_state=random_state)
importance_eval(model, X, y, kfold)

Retornando:

OverallQual: 0.215 +/- 0.013
GrLivArea: 0.152 +/- 0.014
YearBuilt: 0.029 +/- 0.007
BsmtFinSF1: 0.025 +/- 0.004
TotalBsmtSF: 0.023 +/- 0.005
YearRemodAdd: 0.019 +/- 0.003
OverallCond: 0.012 +/- 0.004
GarageCars: 0.012 +/- 0.003
LotArea: 0.009 +/- 0.004
GarageArea: 0.006 +/- 0.003

Podemos notar que, como suspeitamos na análise exploratória dos dados (a qual vamos nos referir a partir de agora por AED), qualidade geral da casa e área útil sobre o nível do solo desempenham papel relevante para determinação do preço da casa. Interessante que informações sobre a garagem (número carros / área) também apresentam relevância para a decisão do GB).

Tratando dados faltantes

Verificamos inicialmente qual o tamanho de nosso problema de missing values por:

def missing_values_check(df):
    missing_val_total = df[df.columns[df.isna().any()].tolist()].isna().sum().sort_values(ascending=False)
    missing_val_perc = missing_val_total/len(df)
    print(pd.concat([missing_val_total, missing_val_perc], axis=1, keys=['Qtde', '%']))

missing_values_check(df)

Que nos indica:

              Qtde         %
PoolQC        1453  0.995205
MiscFeature   1406  0.963014
Alley         1369  0.937671
Fence         1179  0.807534
FireplaceQu    690  0.472603
LotFrontage    259  0.177397
GarageYrBlt     81  0.055479
GarageType      81  0.055479
GarageFinish    81  0.055479
GarageQual      81  0.055479
GarageCond      81  0.055479
BsmtFinType2    38  0.026027
BsmtExposure    38  0.026027
BsmtFinType1    37  0.025342
BsmtCond        37  0.025342
BsmtQual        37  0.025342
MasVnrArea       8  0.005479
MasVnrType       8  0.005479
Electrical       1  0.000685

Assim, começando do mais preocupante para o menos, podemos imaginar que os dados faltantes de PoolQC devem representar inexistência de piscina, já que vimos na AED que o pessoal em Ames não parece gostar muito de piscinas. Vamos checar isso verificando qual a qualidade das piscinas que tem área maior que zero:

print(df[df['PoolArea'] > 0]['PoolQC'].count())

7

Ou seja, as 1453 observações sem dados de qualidade representam casas com área zero de piscina. Assim, criaremos a categoria 'NA' para estas. Faremos o mesmo para MiscFeature que, como podemos verificar pelas categorias em sua distribuição na AED, corresponde a itens especiais diversos (como quadra de tênis), logo é de se esperar que a maioria das casas não os possua. Também o faremos para Fence (provavelmente são casas sem cercas), indicando que tais dados faltantes são na realidade 'NA'.

Já para Alley, provavelmente a falta de valores se deve a casas que não possuem acessos diferenciados (o acesso deve ser pela rua principal, que passa em frente a casa, por exemplo). Assim, também poderemos preencher com valor NA os dados faltantes.

Para FireplaceQu (qualidade das lareiras), provavelmente a falta de dados significa que o número de lareiras é zero. Checamos por:

print(len(df[df['FireplaceQu'].isna()]['Fireplaces'] == 0))

690

Que é exatamente o número de dados faltantes. Então para esses dados faltantes também imputaremos NA.

Para LotFrontage não sabemos se o problema é em função da falta de dados ou inexistência da possibilidade de termos uma medida de frente da propriedade devido a suas características. Como o número de dados faltantes é alto, vamos ignorar esse feature, excluindo-o.

Para GarageYrBlt, GarageType, GarageFinish, GarageQual, GarageCond, notamos que os dados faltantes são em mesma quantidade para todos, o que indica que o problema se deve à falta de garagem nestas casas. Checamos isso cruzando com o área da garagem (igual a zero):

print(len(df[df['GarageType'].isna()]['GarageArea'] == 0))

81

O que confirma nossa hipótese. Logo, também utilizaremos NA para as observações categóricas. Com relação a GarageYrBlt, vamos exclui-la, visto que as outras features de datas parecem possuir melhor poder preditivo.

Assim, imputando NA para todas as features conforme prometido até então, bem como excluíndo as features GarageYrBlt e LotFrontage:

feature_list = ['PoolQC', 'MiscFeature', 'Fence', 'Alley', 'FireplaceQu', 'GarageType', 'GarageFinish', 'GarageQual', 'GarageCond']
for feature in feature_list:
    df[feature] = df[feature].fillna('NA')

df.drop(columns=['GarageYrBlt', 'LotFrontage'], inplace=True)

Os dados faltantes de BsmtFinType2, BsmtExposure, BsmtFinType1, BsmtCond e BsmtQual parecem indicar a inexistência de porão, mas as duas últimas apresentam uma observação a mais sem dados. Assim, para avaliar nossa hipótese, verificando a área de porão disponível:

len(df[df['BsmtFinType1'].isna()]['TotalBsmtSF'] == 0)

37

O que parece confirmar a nossa hipótese de inexistência de porão para as observações em questão. Assim, imputamos NA para essas observações:

feature_list = ['BsmtFinType1', 'BsmtCond', 'BsmtQual']
for feature in feature_list:
    df[feature] = df[feature].fillna('NA')

Como restará um valor valtante para BsmtFinType2 e BsmtExposure, a preenchemos com 'Unf' e 'No', que representam as categorias predominantes nessas features:

df['BsmtFinType2'] = df['BsmtFinType2'].fillna('Unf')
df['BsmtExposure'] = df['BsmtExposure'].fillna('No')

Com relação as 8 observações sem dados de MasVnrArea e MasVnrType: o tipo mais comum para MasVnrType é None, conforme observamos em sua distribuição. Logo, podemos utilizar essa categoria para as 8 observações e lançar como área o valor 0, que é o que ocorre com a maioria das observações.

E temos apenas uma observação com dados faltantes para Electrical, em que preenchemos com a categoria predominante ('SBrkr').

Assim, realizando essas imputações:

df['MasVnrType'] = df['MasVnrType'].fillna('None')
df['MasVnrArea'] = df['MasVnrArea'].fillna(0)
df['Electrical'] = df['Electrical'].fillna('SBrkr')

E por fim verificamos se fomos bem sucedidos (no sentido de não haver mais dados faltantes):

missing_values_check(df)

Que nos retorna:

Empty DataFrame
Columns: [Qtde, %]
Index: []

Indicando que não há mais dados faltantes em nosso dataset.

0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Engenharia de Features e Modelagem (Parte 2)

Criação de Novas Features

Nessa etapa exploramos a criação de algumas novas features, baseadas nas existentes. Por exemplo, possuímos a área útil total acima do solo, a área útil de porão e área da garagem. Podemos criar então uma feature com a área total da casa, a soma dessas três, e podemos verificar que essa nova feature tem uma correlação maior com o preço de venda que as três anteriores:

df['HouseArea'] = df['GrLivArea'] + df['TotalBsmtSF'] + df['GarageArea']
print(df['HouseArea'].corr(df['SalePrice']))

0.8075184760515011

O mesmo faremos para área de varandas: possuímos quatro variáveis de área (OpenPorchSF, EnclosedPorch, 3SsnPorch, ScreenPorch) com pouco ou nenhum poder preditivo, nem correlações significantes com outras variáveis. Então somamos as features em uma nova:

df['PorchArea'] = df['OpenPorchSF'] + df['EnclosedPorch'] + df['3SsnPorch'] + df['ScreenPorch']

Com relação ao número de banheiros, FullBath tem interessante correlação com o preço de venda, como pode ser visto no mapa de calor da AED (0.56). Mas há outras três features relativas a banheiros e lavabos, mas com correlações piores (HalfBath, BsmtFullBath, BsmtHalfBath). Então uma ideia pode ser gerarmos uma nova feature de banheiros, somando todos eles e considerando lavabos como "meio" banheiros. Note que assim criamos uma variável com correlação com o preço de venda maior que todas as medidas individuais:

df['TotalBath'] = df['FullBath'] + 0.5 * df['HalfBath'] + df['BsmtFullBath'] + 0.5 * df['BsmtHalfBath']
print(df['TotalBath'].corr(df['SalePrice']))

0.6317310679319873

Com relação às variáveis de datas, anuais, sabemos que YearRemodAdd (ano da última reforma) tem valor mínimo em 1950, se a casa foi construída antes disso, provavelmente por algum problema com a base de dados. Além disso, sabemos que o ano de venda não parece importar, mas provavelmente importa mais a idade da casa no momento da venda do que a idade da casa em relação a hoje ou a última data disponível na base de dados. Então vamos:

  • Substituir o ano de reforma pelo ano de construção, se a casa foi construída antes de 1950;

  • Criar uma nova feature binária indicando se a casa foi reformada ou não (pela diferença entre ano da última reforma e ano da construção);

  • Criar um feature com o tempo decorrido desde a construção até o momento da venda;

Fazemos isso por:

df.loc[df['YearBuilt'] < 1950, 'YearRemodAdd'] = df[df['YearBuilt'] < 1950]['YearBuilt']
df['IsRemod'] = np.where(df['YearBuilt'] == df['YearRemodAdd'], 0, 1)
df['AgeSold'] =  df['YrSold'] - df['YearBuilt']

Conversão de features númericas para categóricas nominais

Lembramos que MSSubClass, apesar de numérica, é categórica nominal na realidade. Assim a convertemos para nominal. Para a feature do mês de venda não há ordinalidade. Assim também a converteremos para categórica nominal. Faremos o mesmo para o ano de venda:

features_list = ['MSSubClass', 'MoSold', 'YrSold']

for feature in features_list:
    df[feature] = df[feature].astype('category', copy=False)

Agrupamento de categorias em features categóricas nominais

Como vimos na análise exploratória dos dados, o bairro é um dos features categóricos com maior capacidade de diferenciar o preço da casa. Porém, pelo alto número de categorias, e pela distribuição das casas nos bairros, podemos verificar que alguns bairros são pouco representados. Se filtrarmos, por exemplo, os bairros que possuem menos de 2% da observações:

s = df['Neighborhood'].value_counts() < (df.shape[0] * 0.02)
print(s[s].index.tolist())

['ClearCr', 'SWISU', 'StoneBr', 'MeadowV', 'Blmngtn', 'BrDale',
'Veenker', 'NPkVill', 'Blueste']

Vamos "unir" esses bairros nos bairros com distribuição "vizinha" (vamos selecionar aquele com mediana mais próxima), a partir da análise do boxplot da AED, por:

neighborhood_map = [
    ['ClearCr', 'SWISU', 'NPkVill', 'Blueste', 'StoneBr', 'Blmngtn', 'MeadowV', 'BrDale', 'Veenker'],
    ['Crawfor', 'NAmes', 'NAmes',   'NAmes',   'NoRidge', 'CollgCr', 'IDOTRR',  'IDOTRR', 'Somerst']]

for from_, to_ in zip(neighborhood_map[0], neighborhood_map[1]):
    df.loc[df['Neighborhood']==from_, 'Neighborhood'] = to_

Obtendo assim o seguinte boxplot para os bairros apoós o agrupamento:

order = df[['Neighborhood', 'SalePrice']].groupby(by=['Neighborhood']).median().sort_values(by=['SalePrice'], ascending=False).index.tolist()
fig, ax = plt.subplots(figsize=(12, 9))
g = sns.boxplot(x=df['Neighborhood'], y=df['SalePrice'], order=order)
_ = g.set_xticklabels(g.get_xticklabels(), rotation=-25)

A imagem será apresentada aqui.

Raciocínio similar pode ser aplicado para outras features categóricas, em que algumas categorias, de determinada feature, possuem pouca representatividade, e dado algum possível poder preditivo das demais categorias dessa feature, poderia ser interessante proceder com agrupamento. Propomos assim então os seguintes agrupamentos:

  • Foundation: manter as 3 categorias com o maior número de observações (PConc, CBlock e BrkTil) e agrupar as demais;
  • MiscFeature: manter apenas Shed e agrupar as demais;
  • Electrical: manter SBrkr e FuseA e agrupar as demais;
  • Condition1: agrupar PosN, RRNn, PosA, RRNe em 'other' e incluir RRAe com Feedr;
  • HouseStyle: manter 2Story, 1Story, 1.5Fin, Foyer e Slvl e agrupar as demais;
  • GarageType: manter Attchd, Detchd e BuiltiIn e agrupar as demais;
  • SaleType: manter WD, New e COD e agrupar as demais;
  • MSZoning: agrupar RM e RH;
  • Fence: agrupar todas as categorias, com exceção da que representa o valor Good Privacy ("GdPrv");
  • LandSlope: agrupar Mod com Sev;
  • ExterCond: agrupar Ex com Gd e Po com Fa;
  • BsmtCond: agrupar Po com Fa;
  • HeatingQC: agrupar Po com Fa;
  • GarageQual: agrupar Ex com Gd e Po com Fa;
  • GarageCond: agrupar Ex com Gd e Po com Fa;

E implementamos por:

feature_list = ['Foundation', 'MiscFeature', 'Electrical', 'Condition1', 'HouseStyle', 'GarageType', 'SaleType', 'MSZoning', 'Fence', 'LandSlope', 'ExterCond', 'BsmtCond', 'HeatingQC', 'GarageQual', 'GarageCond']

map_ = [
    [['Slab', 'Stone', 'Wood'],
     ['Gar2', 'TenC', 'Othr'],
     ['FuseF', 'FuseP', 'Mix'],
     ['PosN', 'RRNn', 'PosA', 'RRNe'] + ['RRAe'],
     ['1.5Unf', '2.5Unf', '2.5Fin'],
     ['Basment', 'CarPort', '2Types'],
     ['ConLD', 'ConLw', 'ConLI', 'CWD', 'Oth', 'Con'],
     ['RH', 'RM'],
     ['MnPrv', 'GdWo', 'MnWw'],
     ['Sev'],
     ['Ex', 'Po'],
     ['Po'],
     ['Po'],
     ['Ex', 'Po'],
     ['Ex', 'Po']
    ], 
    [['other'] * 3,
     ['NA'] * 3,
     ['other'] * 3,
     ['other'] * 4 + ['Feedr'],
     ['other'] * 3,
     ['other'] * 3,
     ['other'] * 6,
     ['other'] * 2,
     ['other'] * 3,
     ['Mod'],
     ['Gd', 'Fa'],
     ['Fa'],
     ['Fa'],
     ['Gd', 'Fa'],
     ['Gd', 'Fa']
    ]
]

for i, feature in enumerate(feature_list):
    for from_, to_ in zip(map_[0][i], map_[1][i]):
        # print(feature, ":", from_, "->", to_)
        df.loc[df[feature]==from_, feature] = to_

Conversão de features categóricos ordinais para numéricos

Algumas features categóricas são ordenadas, mas seus valores estão no formado de strings. Assim, para aquelas em que temos todos os valores (não há valores como NA, por exemplo), vamos convertê-las para o formato de inteiros que representem a ordenação, conforme concluímos pelos boxplots da análise exploratória:

def encodeOrdinal(feature_list, categories_list):
    for feature, categories in zip(feature_list, categories_list):
        encoder = OrdinalEncoder(categories=[categories])
        encoder.fit(df[feature])
        df[feature] = encoder.fit_transform(df[feature])

feature_list = [['LandSlope'], 
                ['ExterQual'], 
                ['ExterCond'], 
                ['HeatingQC'], 
                ['KitchenQual'], 
                ['PavedDrive']]

slope_cat = ['Gtl', 'Mod']
qual_1_cat = ["NA", "Fa", "TA", "Gd", "Ex"]
qual_2_cat = ["NA", "Po", "Fa", "TA", "Gd", "Ex"]
qual_3_cat = ["NA", "Fa", "TA", "Gd"]
drive_cat = ["N", "P", "Y"]

categories_list = [slope_cat, 
                   qual_1_cat, 
                   qual_3_cat, 
                   qual_1_cat, 
                   qual_2_cat, 
                   drive_cat]

encodeOrdinal(feature_list, categories_list)
0 votos
respondida Dez 15, 2020 por Carlos Alexandre (51 pontos)  
editado Dez 18, 2020 por Carlos Alexandre

Engenharia de Features e Modelagem (Parte 3)

Seleção de Features

Na análise exploratória dos dados verificamos que algumas features possuem um categoria / valor que corresponde a quase totalidade das amostras, não possuindo assim muito poder para diferenciação de preço e tendo alto risco de gerar apenas ruído. Decidimos assim por eliminá-las:

drop_list = ['PoolArea', 'PoolQC', 'Utilities', 'Condition2', 'RoofMatl', 'Heating', 'Street']
df.drop(drop_list, axis=1, inplace=True)

Com criamos algumas features novas, as features utilizadas em sua criação acabam por ser redundantes. Por esse motivo, eliminamos as seguintes features:

drop_list = [
    'TotalBsmtSF', 'GarageArea', 
    'FullBath', 'HalfBath', 'BsmtFullBath', 'BsmtHalfBath',
    'OpenPorchSF', 'EnclosedPorch', '3SsnPorch', 'ScreenPorch',
    'YearBuilt', 'YearRemodAdd'
    ]
df.drop(drop_list, axis=1, inplace=True)

Ainda possuímos algumas variáveis de área causando multicolinearidade perfeita (como verificamos na AED), e assim eliminamos as menos informativas:

drop_list = [
    'BsmtFinSF2', 'LowQualFinSF'
    ]
df.drop(drop_list, axis=1, inplace=True)

E por fim eliminamos aquelas features que não apresentaram correlação, associação ou algum padrão razoável observado nos gráficos da AED relacionado com o preço de venda da casa:

drop_list = [
    'MiscFeature',
    'ExterCond',
    'KitchenAbvGr',
    'MiscVal'
    ]
df.drop(drop_list, axis=1, inplace=True)

Outliers

Um de nossos principais preditores de preço é a área útil da casa. Como observamos no AED, possuímos algumas observações com área acima do solo superior a 4000 pés quadrados que destoam do padrão observado nas demais casas. Assim, optamos por eliminar essas observações:

df = df[df['GrLivArea'] < 4000]

Preparação dos Dados para Modelagem

Para as variáveis categóricas, aplicaremos one-hot-enconding. Após isso, eliminaremos, se existentes, aquelas dummies que possuírem menos que 10 observações.

Para cada feature numérica, verificaremos a medida de assimetria (skewness) e, se maior que 1, indicando considerável assimetria a direita, aplicaremos transformação logarítmica visando reduzi-la. Alguns algoritmos não são sensíveis a tal transformação, mas como estamos utilizando mais de um e podemos incluir facilmente outros nessa implementação, optamos por essa abordagem:

df = pd.get_dummies(df)

X = df.drop(columns = ['Id', 'SalePrice'])
y = np.log(df['SalePrice'])

numercial_features = X.select_dtypes(exclude=['uint8']).columns.tolist()

X_num = X[numercial_features].copy()
X_cat = X.drop(numercial_features, axis=1)

# Tratamento assimetria:
for feature in numercial_features:
    if X_num[feature].skew() > 1:
        # print("before:", df[feature].skew())
        X_num[feature] = np.log(X_num[feature] + 1)
        # print("after:", df[feature].skew())

# Esparsidade:
s = X_cat.sum() < 10
X_cat.drop(columns=s[s].index.tolist(), inplace=True)

# Composição Final:
X = pd.concat([X_num, X_cat], axis=1)

Reavaliação e Seleção de modelos

Após nossa "engenharia de features", testamos novamente os modelos em suas configurações padrão, em validação cruzada (reportamos os scores sobre teste), e é interessante notar como conseguimos melhorar o resultado do modelo de Regressão Linear.

results = evaluate(models, X, y, kfold)

A imagem será apresentada aqui.

Porém, todos os modelos além do de regressão linear podem ser aperfeiçoados pela seleção de seus hiperparâmetros. Vamos avaliar o quanto podemos aperfeiçoá-los em validação cruzada. Para isso, utilizaremos o GridSearchCV e RandomizedSearchCV do sklearn, que buscam, exaustivamente ou aleatoriamente, respectivamente, os hiperparâmetros que maximizam alguma métrica (neg. mean squared error em nosso caso), com base em uma lista de possíveis hiperparâmetros que gostaríamos de avaliar para cada um:

lasso_parameters = {'alpha': np.linspace(0.1, 1.0, num=(10 ** 2)) ** 6}
ridge_parameters = {'alpha': np.linspace(1, 100, num=100)}
gbr_parameters = {
    "loss": ['ls', 'lad', 'huber', 'quantile'],
    "learning_rate": [0.01, 0.02, 0.025, 0.033, 0.05, 0.075, 0.1, 0.125, 0.15, 0.2],
    "min_samples_split": [0.05, 0.075, 0.1, 0.12, 0.15, 0.2, 0.25, 0.5, 1, 2, 5],
    "min_samples_leaf": [0.1, 0.25, 0.5, 1, 2, 5, 10],
    "max_depth": [3, 5, 8],
    "max_features": ["auto", "log2","sqrt"],
    "criterion": ["friedman_mse", "mse"],
    "subsample": [0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0],
    "n_estimators": [250, 500, 1000]
    }

gs_lasso = GridSearchCV(Lasso(), lasso_parameters, cv=kfold, scoring='neg_mean_squared_error')
gs_lasso.fit(X, y)
print(gs_lasso.best_estimator_)

gs_ridge = GridSearchCV(Ridge(), ridge_parameters, cv=kfold, scoring='neg_mean_squared_error')
gs_ridge.fit(X, y)
print(gs_ridge.best_estimator_)

# gs_gbr = GridSearchCV(gbr, parameters, cv=kfold, scoring='neg_mean_squared_error', n_jobs=-1)
gs_gbr = RandomizedSearchCV(GradientBoostingRegressor(random_state=random_state), gbr_parameters, cv=kfold, scoring='neg_mean_squared_error', n_jobs=-1, n_iter=500)
gs_gbr.fit(X, y)
print(gs_gbr.best_estimator_)

E obtemos:

Lasso(alpha=0.00041150149500920354)
Ridge(alpha=19.0)
GradientBoostingRegressor(criterion='mse', learning_rate=0.033, loss='huber',
                          max_depth=5, max_features='sqrt', min_samples_leaf=2,
                          n_estimators=1000, random_state=0, subsample=0.85)

Utilizando o melhor modelo para cada tipo de regressor, realizamos nova validação cruzada:

models = [
    LinearRegression(),
    gs_ridge.best_estimator_,
    gs_lasso.best_estimator_,
    gs_gbr.best_estimator_
    ]

results = evaluate(models, X, y, kfold)

A imagem será apresentada aqui.

E interessante notar que conseguimos resultados, em validação cruzada, muito próximos para os três modelos calibrados, e melhoramos consideravelmente o desempenho do Lasso.

Por fim realizamos um novo exercício de permutation importance para avaliar quais os features com maior impacto no Lasso (para um fold da validação cruzada):

importance_eval(gs_lasso.best_estimator_, X, y, kfold)

O que nos retorna:

HouseArea: 0.134 +/- 0.008
OverallQual: 0.069 +/- 0.007
AgeSold: 0.048 +/- 0.005
GrLivArea: 0.043 +/- 0.004
OverallCond: 0.021 +/- 0.004
LotArea: 0.018 +/- 0.003
TotalBath: 0.005 +/- 0.002

Notar que entre as 7 features mais importantes, 3 foram criadas com base em outras: área total da casa, idade da casa na venda e número total de banheiros. E, como suspeitávamos, a área útil da casa é o fator mais importante na determinação do preço de venda.

Não fizemos a engenharia de features com foco em um estimador em particular, o que poderia ter melhorado o resultado final, em validação cruzada, para algum modelo em específico (que é um exercício também interessante).

Os códigos aqui apresentados estão disponíveis, via notebooks, neste link, e podem ser carregados diretamente no kaggle.

comentou Dez 18, 2020 por Rodrigo Stuckert (66 pontos)  
editado Dez 18, 2020 por Rodrigo Stuckert
Olá Carlos, tudo bem? Gostei bastante da forma que você estruturou seu código, principalmente com a criação de funções que geram gráficos comparando as performances dos diversos modelos testados; inclusive pretendo implementar algumas dessas funções em meu projeto também. Também achei interessante alguns pacotes que você usou, como o dython, que eu desconhecia até então.

Gostaria apenas de fazer algumas considerações quanto a alguns eventuais pontos que provavelmente passaram batidos na hora de transcrever o código pro Prorum, o que é comum, ainda mais tratando-se de um projeto mais extenso. Em alguns pedaços do código, aconteceu alguma coisa que fez com que "ax = ax[i][1]" se tornasse "ax = ax[i][3]", ou outros valores no último dígito. Também, mais ao final do código, mais especificamente na Engenharia de Features parte 3, o pedaço de código "results = evaluate(models, X, y, kfold)" retorna a mensagem "posx and posy should be finite values" e gera gráficos vazios. Acredito que deva ter sido algum erro de digitação ocorrido. No mais excelente trabalho, e parabéns pelo projeto! Forte abraço!
comentou Dez 18, 2020 por Carlos Alexandre (51 pontos)  
Obrigado pela suas observações Rodrigo! Já corrigi os ax's, e com relação ao último erro, é em função de uma linha de código que faltou na cópia aqui no prorum, no segmento de código de "Preparação dos Dados para Modelagem".  Novamente, muito obrigado pelas considerações! Abs.
...