<jupyter_start><jupyter_text>Summarise the dataset<jupyter_code>arabica_df=pd.read_csv('/Users/peace/standalone_projects/experian_coffee/data/arabica_final.csv')
arabica_df=arabica_df.rename(columns={'Country.of.Origin':'Country_of_Origin','Farm.Name':'Farm_Name','Lot.Number':'Lot_Number','ICO.Number':'ICO_Number','Number.of.Bags':'Number_of_Bags','Bag.Weight':'Bag_Weight','In.Country.Partner':'In_Country_Partner','Harvest.Year':'Harvest_Year','Grading.Date':'Grading_Date','Processing.Method':'Processing_Method','Clean.Cup':'Clean_Cup','Total.Cup.Points':'Total_Cup_Points','Category.One.Defects':'Category_One_Defects','Category.Two.Defects':'Category_Two_Defects','Certification.Body':'Certification_Body','Certification.Address':'Certification_Address','Certification.Contact':'Certification_Contact','Cupper.Points':'Cupper_Points'})
arabica_df.describe(include='object')
arabica_df.describe()<jupyter_output><empty_output><jupyter_text>* The minimum values of some cols (Aroma, Flavor, Aftertaste...) are zero, which is quite strange (del the row)<jupyter_code># Delete the only row that has all zero gradings
error_index=arabica_df[arabica_df['Flavor']==0].index.values[0]
arabica_df=arabica_df.drop(error_index)<jupyter_output><empty_output><jupyter_text>Check Missing Values<jupyter_code>total=arabica_df.isnull().sum().sort_values(ascending=False)
percentage=(total/arabica_df.isnull().count()).sort_values(ascending=False)
percentage.head(17)<jupyter_output><empty_output><jupyter_text>Impute Missing Values  QuakersOne missing value* Only 1/54 coffee beans from the same owner have one Quaker<jupyter_code>arabica_df['Quakers']=arabica_df['Quakers'].fillna(0)<jupyter_output><empty_output><jupyter_text>Cupper PointsThe newly scrapped data dosen't have cupper_points col* Fill nan with the median value from the same In_Country_Partner* Fill the rest nan with the median value of the dataset<jupyter_code>arabica_df['Cupper_Points']=arabica_df.groupby('In_Country_Partner').Cupper_Points.apply(lambda s: s.fillna(s.median()))
arabica_df['Cupper_Points']=arabica_df['Cupper_Points'].fillna(arabica_df['Cupper_Points'].median())<jupyter_output><empty_output><jupyter_text>Clustering* Delete Altitude col (>15% missing values) and Number of bags col (irrelavent)* K-means clustering, K-means clustering with cosine similarities, DBSCAN<jupyter_code># Prepare the data
X = arabica_df.select_dtypes(exclude=['object'])
X.drop(['Altitude','Number_of_Bags'],1,inplace=True)
x_columns=X.columns
X.isnull().sum()<jupyter_output><empty_output><jupyter_text>K-Means with euclidean distance<jupyter_code># Elbow method to determine the number of clusters 
from sklearn.cluster import KMeans

clusters=[]

for i in range(1,11):
    km=KMeans(n_clusters=i).fit(X)
    clusters.append(km.inertia_)

fig,ax=plt.subplots(figsize=(7,4))
sns.lineplot(x=list(range(1,11)),y=clusters,ax=ax)
ax.set_title('Inertia of Euclidean k-Means versus number of clusters')
ax.set_xlabel('Clusters')
ax.set_ylabel('Inertia')

plt.show()<jupyter_output><empty_output><jupyter_text>As can be seen from the graph, 3 is the last bend.  K-means with cosine similarities<jupyter_code>from sklearn import preprocessing

clusters=[]

normalized_vectors=preprocessing.normalize(X)
for i in range(1,11):
    km=KMeans(n_clusters=i).fit(normalized_vectors)
    clusters.append(km.inertia_)

fig,ax=plt.subplots(figsize=(7,4))
sns.lineplot(np.arange(1,11),clusters)


plt.xlabel('Number of clusters')
plt.ylabel("Inertia")
plt.title("Inertia of Cosine k-Means versus number of clusters")

plt.show()<jupyter_output><empty_output><jupyter_text>DBSCAN (based on the density of data points)* minPts: at least the number of features in your data* eps:<jupyter_code># Find the optimum eps at the last bend in the graph
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors

neigh = NearestNeighbors(n_neighbors=2)
nbrs = neigh.fit(X)
distances, indices = nbrs.kneighbors(X)
distances = np.sort(distances, axis=0)
distances = distances[:,1]
fig,ax=plt.subplots(figsize=(7,4))
plt.plot(distances)<jupyter_output><empty_output><jupyter_text>Find the eps=3  Evaluation (Sihouette Score)(b-a)/max(a,b)* a- mean intra-cluster distance* b- mean nearest-cluster distance<jupyter_code>from sklearn.metrics import silhouette_score

# prepare models
kmeans=KMeans(n_clusters=3).fit(X)

normalized_vectors=preprocessing.normalize(X)
normalized_kmeans=KMeans(n_clusters=3).fit(normalized_vectors)

min_samples=X.shape[1]+1
dbscan=DBSCAN(eps=3,min_samples=min_samples).fit(X)

print('kmeans: {}'.format(silhouette_score(X, kmeans.labels_, 
                                           metric='euclidean')))
print('Cosine kmeans:{}'.format(silhouette_score(normalized_vectors,
                                          normalized_kmeans.labels_,
                                          metric='cosine')))
print('DBSCAN: {}'.format(silhouette_score(X, dbscan.labels_, 
                                           metric='cosine')))<jupyter_output>kmeans: 0.47903527156227416
Cosine kmeans:0.7617933045417482
DBSCAN: 0.8259323420593372<jupyter_text>As can be seen above, DBSCAN > Cosine > kmeans. Cosine similarity works better due to the high dimensionality of the data  Visualizing Clusters (PCA & t_SNE)Use dimensionality reduction to visualise the data in 2d and 3d<jupyter_code>from sklearn.decomposition import PCA
def pca_prepare(n_components,df,labels):
    pca=PCA(n_components=n_components)
    pca_data=pca.fit_transform(df)
    pca_data=np.vstack([pca_data.T,labels]).T
    pca_df=pd.DataFrame(pca_data,columns=['x','y','labels'])
    return pca_df
pca_df_kmeans=pca_prepare(2,X,kmeans.labels_)
pca_df_cosine_kmeans=pca_prepare(2,X,normalized_kmeans.labels_)
pca_df_DBSCAN=pca_prepare(2,X,dbscan.labels_)

fig,axes=plt.subplots(3,figsize=(10,6))
sns.scatterplot(ax=axes[0],x=pca_df_kmeans.x, y=pca_df_kmeans.y,hue=pca_df_kmeans.labels, palette='Set2').set(title='Kmeans with euclidean')
sns.scatterplot(ax=axes[1],x=pca_df_cosine_kmeans.x, y=pca_df_cosine_kmeans.y,hue=pca_df_cosine_kmeans.labels, palette='Set2').set(title='Kmeans with cosine')
sns.scatterplot(ax=axes[2],x=pca_df_DBSCAN.x, y=pca_df_DBSCAN.y,hue=pca_df_DBSCAN.labels, palette='Set2').set(title='DBSCAN')

fig.suptitle('PCA visualization of the clustering results')
fig.show()<jupyter_output><empty_output><jupyter_text>Only one cluster is identified with DBSCAN methods; -1 label is the label for noises<jupyter_code>from sklearn.manifold import TSNE
def tsne_prepare(n_components,df):
    tsne=TSNE(n_components=n_components)
    tsne_data=tsne.fit_transform(df)
    return tsne_data
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import ListedColormap
tsne_data=tsne_prepare(3, X)

tsne_data_kmeans=np.vstack([tsne_data.T,kmeans.labels_]).T
tsne_3d_df_kmeans=pd.DataFrame(tsne_data_kmeans,columns=['x','y','z','labels'])

tsne_data_cosine_kmeans=np.vstack([tsne_data.T,normalized_kmeans.labels_]).T
tsne_3d_df_cosine_kmeans=pd.DataFrame(tsne_data_cosine_kmeans,columns=['x','y','z','labels'])

tsne_data_dbscan=np.vstack([tsne_data.T,dbscan.labels_]).T
tsne_3d_df_dbscan=pd.DataFrame(tsne_data_dbscan,columns=['x','y','z','labels'])

fig=plt.figure(figsize=(6,6))
ax=Axes3D(fig)

sc=ax.scatter(tsne_3d_df_cosine_kmeans.x, tsne_3d_df_cosine_kmeans.y,tsne_3d_df_cosine_kmeans.z,c=tsne_3d_df_cosine_kmeans.labels,marker='o')<jupyter_output><empty_output><jupyter_text>Interpreting ClustersFind out which features are more significant in the clustering process<jupyter_code>from sklearn.preprocessing import MinMaxScaler

scaler=MinMaxScaler()
df_scaled=pd.DataFrame(scaler.fit_transform(normalized_vectors),columns=x_columns)

df_scaled['kmeans_cosine']=normalized_kmeans.labels_

df_mean=(df_scaled.groupby('kmeans_cosine').mean())
df_mean
results = pd.DataFrame(columns=['Variable', 'Var'])
for column in df_mean.columns[0:]:
    results.loc[len(results), :] = [column, np.var(df_mean[column])]
selected_columns = list(results.sort_values('Var', ascending=False).head(10).Variable.values) + ['kmeans_cosine']
tidy = df_scaled[selected_columns].melt(id_vars='kmeans_cosine')
plt.figure(figsize=(10,5))
sns.barplot(x='kmeans_cosine', y='value', hue='variable', data=tidy)
plt.title('Features that distinguish clusters')
plt.show()<jupyter_output><empty_output><jupyter_text>Random Forest Feature SelectionUse classification to understand the importance of each feature (More work could be done here for validation and finding the most suitable model.)<jupyter_code># Box plot to identify the trend
arabica_df['group']=normalized_kmeans.labels_
fig,axes=plt.subplots(3,figsize=(8,8))
sns.boxplot(x='group',y='Aroma',data=arabica_df,ax=axes[0])
sns.boxplot(x='group',y='Flavor',data=arabica_df,ax=axes[1])
sns.boxplot(x='group',y='Sweetness',data=arabica_df,ax=axes[2])
# Get Correlation Matric (Heatmap)
corrmat=X.corr()
plt.figure(figsize=(6,6))
sns.heatmap(corrmat)<jupyter_output><empty_output><jupyter_text>Aroma, flavor, aftertaste, acidity, body and balance are positively correlated; uniformity, clean_cup and sweetness are mostly 10/10. Total_Cup_Points are negatively correlated with category_one_defects and category_two defects as expected.<jupyter_code>k=5
cols=corrmat.nlargest(k,'Total_Cup_Points')['Total_Cup_Points'].index
cm=np.corrcoef(X[cols].values.T)
hm=sns.heatmap(cm,cbar=True,annot=True,square=True,fmt='.2f',yticklabels=cols.values,xticklabels=cols.values)<jupyter_output><empty_output><jupyter_text>Flavor and aftertaste are the most positively correlated features to Total_Cup_Points<jupyter_code># Choosing the ideal hyper-parameters (Random SearchCV)
from sklearn.model_selection import RandomizedSearchCV
from sklearn.ensemble import RandomForestClassifier

y=normalized_kmeans.labels_

n_estimators = [int(x) for x in np.linspace(start = 80, stop = 1000, num = 10)]
max_features = ['auto', 'sqrt']
max_depth = [int(x) for x in np.linspace(10, 110, num = 11)]
max_depth.append(None)
min_samples_split = [2, 5, 10]
min_samples_leaf = [1, 2, 4]
bootstrap = [True, False]

random_grid = {'n_estimators': n_estimators,
               'max_features': max_features,
               'max_depth': max_depth,
               'min_samples_split': min_samples_split,
               'min_samples_leaf': min_samples_leaf,
               'bootstrap': bootstrap}

random_forest=RandomForestClassifier()
random_forest_search=RandomizedSearchCV(estimator=random_forest,param_distributions=random_grid,n_iter=100,cv=5,verbose=2,random_state=42,n_jobs=-1)
random_forest_search.fit(X,y)
# Finding the best params
random_forest_search.best_params_
# Evaluate random search
from sklearn.model_selection import KFold, cross_val_score
n_folds=5
def cross_v(model):
    kf=KFold(n_folds,shuffle=True,random_state=42).get_n_splits(X)
    score=cross_val_score(model,X,y,cv=kf)
    return score.mean()
    
base_model=RandomForestClassifier(random_state=42,n_estimators=100)
rf=random_forest_search.best_estimator_

base_score=cross_v(base_model)
rf_score=cross_v(rf)

print(f'Base Model Score: {base_score}')
print(f'Tuned Model Score: {rf_score}')
# Choosing the ideal hyper-parameters (GridSearchCV)
from sklearn.model_selection import GridSearchCV
param_grid = {
    'bootstrap': [False],
    'max_depth': [60, 70, 80, 90],
    'min_samples_leaf': [1, 2, 3],
    'max_features': ['sqrt'],
    'min_samples_split': [3, 5, 7,9],
    'n_estimators': [400, 600,800, 1000]
}

random_forest_1=RandomForestClassifier()
grid_search=GridSearchCV(estimator=random_forest_1,param_grid=param_grid,cv=5,n_jobs=-1,verbose=2)
grid_search.fit(X,y)
# Finding the best params
grid_search.best_params_
# Evaluate grid search
rf_grid=grid_search.best_estimator_

rf_score=cross_v(rf)
rf_grid_score=cross_v(rf_grid)

print(f'Random Tuned Score: {rf_score}')
print(f'Grid Tuned Score: {rf_grid_score}')<jupyter_output>Random Tuned Score: 0.9665819783197833
Grid Tuned Score: 0.96937185443283<jupyter_text>It seems we have about maxed out performance.<jupyter_code># Feature Importance Visualization

data=np.array([rf.feature_importances_,x_columns]).T
columns=list(pd.DataFrame(data,columns=['Importance','Feature']).sort_values('Importance',ascending=False).head(10).Feature.values)
tidy=df_scaled[columns+['kmeans_cosine']].melt(id_vars='kmeans_cosine')
plt.figure(figsize=(10,5))
sns.barplot(x='kmeans_cosine',y='value',hue='variable',data=tidy)
plt.title('Features that distinguish clusters')
plt.show()<jupyter_output><empty_output><jupyter_text>* Label 2 coffee beans seem to have the highest number of defects and the lowest grading in terms of aftertaste, flavor, and aroma, which could possibly suit James who prefers coffee the american way (if low gradings in flavor and aroma mean less flavor and aroma).* In regards to the sweetness, coffee beans from all three cluster groups have very high gradings.To find the best coffee bags that would satisfy everyone in the team, I'd look into the cluster group 1 and 2, and find the sweetest coffee beans with the least defects made in Brazil, preferrably from South of Minas. As the only coffee beans belong to group 2 and was processed in Brazil has 55 category two defects, I choose to mainly consider the ones from group 1.<jupyter_code>final_group=arabica_df[(arabica_df['Group']==1)&(arabica_df['Country_of_Origin']=='brazil')&(arabica_df['Sweetness']==10)&(arabica_df['Category_One_Defects']==0)]
final_brazil=final_group.sort_values(['Flavor','Aroma','Category_Two_Defects'],ascending=['True','True','True']).head(3)
final_south_of_minas=final_group[final_group['Region']=='south of minas'].sort_values(['Flavor','Aroma','Category_Two_Defects'],ascending=['True','True','True']).head(3)<jupyter_output><empty_output><jupyter_text>The below coffee beans could potentially be the ones that help Marcio achieve his ambitious goals! The first 3 are the ones that have the lowest grading of aroma, flavor and the lowest number of category two defects, and the last three are the ones from south of Minas with similar gradings.<jupyter_code>final_brazil
final_south_of_minas<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1258855199773307, 'fraction_numerical': 0.018560498724851233, 'mean_word_length': 8.05516356638871, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 78, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '422774', 'license': 'MIT', 'n_tokens_mistral': 5373, 'n_tokens_neox': 5041, 'n_words': 885, 'path': 'Exploratory Analysis.ipynb', 'repo_name': 'Phil-avi/coffee-quality-database', 'size': '14116'}