Os dados disponibilizados no site são referentes a 137 restaurantes distintos na base de treino e 100 mil restaurantes na base de teste. A target do problema é receita e, como não temos dados sobre receita na base de teste, iremos trabalhar somente com a base de treino
O primeiro problema que temos ao analisar os dados é o tamanho da base de treino. 137 observações e 43 colunas é receita para overfit, pois as poucas observações disponíveis não permitem que o modelo seja generalizado. Em outras palavras, fica muito mais fácil encaixar modelos com uma amostragem de treino pequena.
O segundo ponto que chama atenção é que 37 das variavés disponíveis são anônimas. Não sabemos o que elas representam. Apenas que dentre elas temos: população em determinada area, idade, distribuição de gênero, custo do metro quadrado, disponibilidade de estacionamento, presença de escolas, bancos. São, em geral, dados demográficos, geográficos e comerciais. Já podemos esperar alta multicolinearidade nessa base.
No primeiro momento, após importar as bibliotecas e os dados, vamos tentar conhecer melhor nossa base. Temos variáveis categóricas (essas especificadas) e contínuas (anônimas). Vamos começar pelas categóricas. São 4 no total, uma referente ao ano de inauguração da loja, uma indicando a cidade onde está localizada, uma indicando se está em região de metropole ou não e a última indica o tipo de restaurante.
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
dataset = pd.read_csv('/home/gusta/Área de Trabalho/Restaurant pred compt/train.csv')
dataset.dtypes
dataset.isnull().sum()
dataset = dataset.set_index('Id')
sns.factorplot(x='City Group',data=dataset,kind='count',size=3,aspect=0.75)
plt.title('Quantidade de cidades')
dataset[["City Group","revenue"]].groupby(["City Group"]).mean().plot(kind="bar")
plt.title('Média de receita')
dataset[["City","revenue"]].groupby(["City"]).mean().plot(kind="bar")
plt.title('Receita Media Cidades')
sns.countplot(x='Type', data=dataset)
plt.title('Distribuição Tipos')
O que podemos inferir destar variáveis categóricas? Fica claro que restaurantes situados em metrópoles tem receita maior. Vamos transformar essa primeira coluna em dummie. Agora já em relação às cidades, não podemos inferir nada. Na verdade, já tendo a informação sobre metropole e as possíveis variáveis geograficas e demograficas escondidas nas outras colunas não faz sentido manter a coluna de cidade.
Outra coluna categórica que chamou atenção foi a que define o tipo de restaurante. No total existem 4 tipos de restaurante possíveis, porem na amostra dedicada ao treino temos a clara dominâcia de duas delas, enquanto uma categoria tem apenas uma observação e a outra restante não tem nenhuma. Vamos eliminar essa coluna pois não queremos manter uma variável desbalanceada para nosso treino. A última coluna com informação categórica diz respeito à data de inauguração do restaurante. Vamos usar essa informação, mas não do jeito que está disponível. Iremos utilizar a data que os dados foram disponibilizados como o dia em que as amostras foram colhidas e iremos, a partir desta data, calcular o tempo, em anos, que o restaurante está aberto. Segue o código
dataset = pd.get_dummies(dataset, columns = ['City Group'])
dataset = dataset.rename(columns={'City Group_Big Cities':'Big City', 'City Group_Other':'Other'})
dataset['Open Date'] = pd.to_datetime(dataset['Open Date']).dt.year
dataset['Hoje'] = pd.DataFrame({'Date':np.repeat(['01/01/2015'],[len(dataset)]) })
dataset['Hoje'] = pd.to_datetime(dataset['Hoje']).dt.year
dataset['Anos Aberto'] = dataset['Hoje'] - dataset['Open Date']
dataset = dataset.drop(['City', 'Type', 'Other', 'Open Date', 'Hoje'], axis = 1 )
Agora temos somente variáveis contínuas e uma dummie. Vamos analisar o mapa de correlação como todo. Isso pode nos ajudar a entender as variáveis anônimas. Vamos tambem buscar informações mais claras da relação direta com o \( \it Target\) e a distribuição algumas colunas anônimas (apesar do código gerar infomação sobre a distribuição de todas colunas anônimas, achei melhor não colocar aqui para não poluir com muitas imagens).
plt.figure(figsize=(50,30))
corr = dataset.corr()
sns.heatmap(corr,
xticklabels=corr.columns.values)
corr_revenue = dataset.corr()['revenue'].sort_values(ascending=False)
plt.figure(figsize=(10,7))
corr_revenue.drop('revenue').plot.bar()
plt.show()
for i in range(1,38):
plt.subplots(1, figsize=(5, 2))
sns.distplot(dataset['P{}'.format(i)], kde=False);
plt.show()
Como era de se esperar, existe grande correlação entre várias colunas. Uma forma de reduzirmos a quantidade de colunas e ao mesmo tempo trabalhar para diminuir a correlação entre elas, reduzindo a possibilidade de overfitting, é por meio do \(\it Principal\) \(\it Component\) \(\it Analysis\), o PCA.
Tal metódo consiste, de forma resumida, em transformar um conjunto de variáveis correlacionadas em um conjunto menor, chamados de componentes principais. Os componentes principais são projeções ortogonais das variáveis originais e são ordenados de acordo com sua contribuição para explicar a variação nos dados. E é essa variância que é usada para informação, ou seja, para definir em quantos componentes as colunas serão transformadas.
Como se trata de dados diversos, trazendo informações que com certeza não estão no mesmo padrão, iremos padronizar os features para evitar que o PCA seja prejudicado caso haja variâncias muito grandes ou pequenas. Usaremos a ferramenta StandartScale do sklearn para esse processo. Note que vamos usar apenas os dados do treino para moldar os parâmetros da uniformização, mantendo os dados de teste isolados. Para isso vamos separar a base entre \(\it features\) e \(\it target\) e usar 15% dos dados para teste. Note que não estamos normalizando a última coluna, que é a única dummie da base.
X = dataset.iloc[:, np.r_[0:37, 38:40]].values
y = dataset.iloc[:, 37].values
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.15)
SCx = StandardScaler()
X_train[:, :-1] = SCx.fit_transform(X_train[:, :-1])
X_test[:, :-1] = SCx.transform(X_test[:, :-1])
Agora vamos ao PCA. Novamente iremos recorrer ao sklearn para os próximos passos. No primeiro momento vamos fazer a função do PCA sem parâmetros, apenas iremos usar a base de treino das features e iremos calcular o número de componentes para explicar a variância. Vamos mostrar a contribuição de cada componente e mostrar o acumulado da variancia explicada também. Nesse momento estabelecemos até onde queremos a variância acumulada e escolhemos o número de componentes baseado nessa informação.
from sklearn.decomposition import PCA
pca = PCA()
pca.fit(X_train) #Apenas fit. Transform vem depois de achar o número de componentes
cumulative_explained_variance = np.cumsum(pca.explained_variance_ratio_)
plt.plot(range(len(cumulative_explained_variance)), cumulative_explained_variance)
plt.bar(range(1, (len(pca.explained_variance_ratio_) +1)), pca.explained_variance_ratio_, alpha=0.5, align='center',
label='Explicação individual')
plt.step(range(1, (len(pca.explained_variance_ratio_) +1)), np.cumsum(pca.explained_variance_ratio_), where='mid',
label='Variância explicada acumulada')
plt.ylabel('Taxa de Variância Explicada')
plt.xlabel('Nº Componentes')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
variance = pd.DataFrame(pca.explained_variance_ratio_)
np.cumsum(pca.explained_variance_ratio_)
Após rodas o código acima, concluímos que o número de componentes para explicar a variação de 90 da variância é 9. Agora vamos treinar novamente a função PCA , pórem ocorrerá a transformação dos dados também. Vale frisar que o \( \it fit\) ocorre somente na base de treino, tanto na padronização quanto na aplicação do PCA. Isso evita que ocorra contaminação da base de teste. Ambos os processos também só ocorrem nos \(\it features\). A \(\it Target\) não é tocada em nenhum dos momentos.
pca = PCA(n_components = 9)
X_train = pca.fit_transform(X_train)
X_test = pca.transform(X_test)
Agora com nossos \( \it features\) transformados, vamos verificar qual modelo é mais adequado.
Para isso vamos utilizar a técnica de Cross Validation. CV consiste em dividir a base de treino em k subgrupos de mesma dimensão. Após a divisão, usamos k-1 subgrupos para treinar o modelo e testar no subgrupo que foi deixado de fora. Tal processo é feito repetivamente até todos k subgrupos tenham sidos usados como teste pelo menos uma vez, ou seja, a quantidade de \( \it output\) vai depender de quantos subgrupos foram formados. Nesse processo de validação pode haver um otimismo nos resultados. Isso porque, diferentemente da base de teste separada no começo do código, que não teve nenhum contato com os processos elaborados na base de treino, os subgrupos usados no CV foram "contaminados" pela base de treino.
from sklearn.model_selection import KFold, cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
kfold = KFold(n_splits=5)
kfold.get_n_splits(X_train)
print(kfold)
for train_index, test_index in kfold.split(X_train):
X_train_val, X_test_val = X[train_index], X[test_index]
y_train_val, y_test_val = y[train_index], y[test_index]
mse_medias, mse_mins, mse_maxs = list(),list(),list()
r2_medias, r2_mins, r2_maxs = list(),list(),list()
models = {LinearRegression(), DecisionTreeRegressor()}
def CV(model):
scores = cross_val_score(model, X_train_val, y_train_val, scoring='r2', cv=kfold)
mse = cross_val_score(model, X_train_val, y_train_val, scoring='neg_mean_squared_error', cv=kfold)
r_mean, r_min, r_max = scores.mean(), scores.min(), scores.max()
m_mean, m_min, m_max = mse.mean(), mse.min(), mse.max()
r2_medias.append(r_mean)
r2_mins.append(r_mean - r_min)
r2_maxs.append(r_max - r_mean)
mse_medias.append(m_mean)
mse_mins.append(m_mean - m_min)
mse_maxs.append(m_max - m_mean)
for model in models:
CV(model)
Modelos = ['Linear Reg','DT', 'KNN', 'RF' ]
CV_r2 = pd.DataFrame(list(zip(Modelos, r2_mins, r2_medias, r2_maxs)),
columns=['Modelo','R2 Min', 'R2 Media', 'R2 Max'])
CV_mse = pd.DataFrame(list(zip(Modelos, mse_mins, mse_medias, mse_maxs)),
columns=['Modelo','MSE Min', 'MSE Media', 'MSE Max'])
sns.factorplot(y='Modelo',x='R2 Max',data=CV_r2,kind='bar',size=5,aspect=2)
sns.factorplot(y='Modelo',x='MSE Min',data=CV_mse,kind='bar',size=5,aspect=2)
Após todo o processo, definimos o melhor modelo baseado no R2 máximo e menor Mean Squared Error. O modelo escolhido é o KNN.
def KNN():
RegKNN = KNeighborsRegressor()
RegKNN.fit(X_test, y_test)
y_pred = RegKNN.predict(X_test)
plt.scatter(X_test[:,0], y_test, color = 'red')
plt.scatter (X_test[:,0], y_pred, color = 'blue')
plt.ylabel('Receita')
plt.xlabel('Restaurantes')
plt.title('KNN - n = 5')
plt.show()
print (r2_score( y_test, y_pred))
KNN()
Vemos que a previsão está longe de ser satisfatória (score = 0.086). Era esperado dado o resultado no CV, apesar do maior entre os testados. A falta de disponibilidade de dados também prejudicou o modelo. Uma base de treino maior e mais heterogênea poderia contribuir. Além disso, informações claras sobre quais informações as colunas trazem facilitariam o processo de preparo e escolha dos \( \it features\).