{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy.stats import ttest_ind\nfrom sklearn import metrics\nfrom scipy import stats\n\n\n#Importing libraries for model creation\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.mixture import BayesianGaussianMixture\n\n#Importing pre-processing\nfrom sklearn import preprocessing\n#Decomposition\nfrom sklearn.decomposition import PCA\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-14T12:17:21.724670Z","iopub.execute_input":"2022-07-14T12:17:21.725129Z","iopub.status.idle":"2022-07-14T12:17:21.739963Z","shell.execute_reply.started":"2022-07-14T12:17:21.725095Z","shell.execute_reply":"2022-07-14T12:17:21.737843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Problem Approach\n\nIn this notebook we're tasked with clustering an unlabeled dataset with the evaluation metric being the RandScore. For our model we compare a Gaussian Mixture with a Bayesian Gaussian mixture and end up setting on the Bayesian Gaussian for the final submission. The overall approach is as follows:\n\n* Data Loading\n* Data Normalization\n* Feature Selection\n* Determining Optimal Number of Clusters\n* Training Gaussian Mix model and Bayesian Gaussian Mix and compare results\n* Submitting result from best model","metadata":{}},{"cell_type":"markdown","source":"## Loading Data","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('../input/tabular-playground-series-jul-2022/data.csv', index_col=False)\ndf = df.fillna(0)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:17:21.778601Z","iopub.execute_input":"2022-07-14T12:17:21.779055Z","iopub.status.idle":"2022-07-14T12:17:22.801230Z","shell.execute_reply.started":"2022-07-14T12:17:21.779017Z","shell.execute_reply":"2022-07-14T12:17:22.799940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Scaling\n\nA Gaussian Mixture model assumes that each variable follows a gaussian distribution. Examining the dataset below we can see this isnt the case with many variables being skewed. In order to correct this I used the  sklearn power transformer implementing the yeo-johnson method, this method doesnt require all datapoints to be positive (like a Box-Cox transform) and resulted in a better score than other transformerms.","metadata":{}},{"cell_type":"code","source":"#Shaping to appropriate format\ndf_copy = df.drop(columns = ['id'])\n\nscaler = preprocessing.PowerTransformer(method = 'yeo-johnson', standardize=True).fit(df_copy.values)\nscaled_df = pd.DataFrame(scaler.transform(df_copy.values), index = df_copy.index, columns = df_copy.columns)\np_vals = []\nfor col in df_copy.columns:\n    pre_transform = stats.shapiro(df[col]).pvalue\n    post_transform = stats.shapiro(scaled_df[col]).pvalue\n    p_vals.append([col, pre_transform, post_transform])\n\np_val_df = pd.DataFrame(p_vals, columns = ['Variable', 'Pre-Transform', 'Post-Transform'])\nprint(p_val_df.sort_values(by=['Pre-Transform']))\n\n\nmelted_df_pre = df_copy.melt(value_vars = df_copy.columns,\n                    value_name = 'Value', var_name = 'Variable')\nmelted_df_post = scaled_df.melt(value_vars = df_copy.columns,\n                    value_name = 'Value', var_name = 'Variable')\nmelted_df_pre['Transform'] = 'No Transform'\nmelted_df_post['Transform'] = 'yeo-johnson'\nmelted_df = pd.concat([melted_df_pre, melted_df_post], ignore_index = True)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:17:22.803474Z","iopub.execute_input":"2022-07-14T12:17:22.803918Z","iopub.status.idle":"2022-07-14T12:17:28.105354Z","shell.execute_reply.started":"2022-07-14T12:17:22.803881Z","shell.execute_reply":"2022-07-14T12:17:28.103968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#melted_df.head(n = 10)\nsns.set(rc = {'figure.figsize':(15,12)})\nv = sns.FacetGrid(melted_df, col='Variable', hue = 'Transform', height=2.5, col_wrap=5, sharex = False)\nv.map(sns.histplot, 'Value', alpha = 0.5).add_legend()\nv.tight_layout","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:17:28.107507Z","iopub.execute_input":"2022-07-14T12:17:28.107913Z","iopub.status.idle":"2022-07-14T12:18:46.208831Z","shell.execute_reply.started":"2022-07-14T12:17:28.107878Z","shell.execute_reply":"2022-07-14T12:18:46.207523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Filtering for Best Variables\n\nBasing the approach off the notebook here: https://www.kaggle.com/code/ricopue/tps-jul22-clusters-and-lgb\n\nwe're going to take a subset of our factors to use for model training.\n\n\n","metadata":{}},{"cell_type":"code","source":"sns.set(rc = {'figure.figsize':(15,12)})\nsns.set_style('white')\nheatmap = sns.heatmap(scaled_df.corr(), annot=False, cmap='BrBG',)\nheatmap.set_title('Variable Correlation', fontdict={'fontsize':26}, pad=16);\n","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:18:46.211606Z","iopub.execute_input":"2022-07-14T12:18:46.212046Z","iopub.status.idle":"2022-07-14T12:18:47.309054Z","shell.execute_reply.started":"2022-07-14T12:18:46.212010Z","shell.execute_reply":"2022-07-14T12:18:47.307762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_data =['f_07','f_08', 'f_09', 'f_10','f_11', 'f_12', 'f_13', 'f_22','f_23', 'f_24', 'f_25','f_26','f_27', 'f_28']\nscaled_df = scaled_df[best_data]","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:18:47.310946Z","iopub.execute_input":"2022-07-14T12:18:47.311804Z","iopub.status.idle":"2022-07-14T12:18:47.322494Z","shell.execute_reply.started":"2022-07-14T12:18:47.311762Z","shell.execute_reply":"2022-07-14T12:18:47.320971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Parameter Tuning\n\nThe first parameter we need to identify is the number of clusters to predict. We'll do this by taking a subset of the data(to reduce training time) and training a model for between 2-15 clusters. \n\nTo compare model performance we'll use the Silhouette score explained here: https://en.wikipedia.org/wiki/Silhouette_(clustering). \n\nGiven that increasing the number of groups will naturally lead to a lower silhouette score we'll use the elbow method explained here: https://en.wikipedia.org/wiki/Elbow_method_(clustering) to look at when the rate of change reduces as we increase the number of clusters.\n\nThe best leaderboard score was a result of using n = 7 clusters","metadata":{}},{"cell_type":"code","source":"sample_df = scaled_df.sample(n = 5000)\nclusters = range(2,15)\nscores = []\n\nfor i in clusters:\n    gm = GaussianMixture(n_components=i, n_init=5, init_params='kmeans',\n                        verbose = 0)\n    gm_prediction = gm.fit_predict(sample_df)\n    # Calculate Silhoutte Score and append to a list\n    score = metrics.silhouette_score(sample_df, gm_prediction, metric='euclidean')\n    scores.append(score)\n    print('Number of Clusters: ', i, ' Score: ', score)\n  \n\nplt.plot(clusters, scores, 'bo-')\nplt.xlabel('Number of Clusters')\nplt.ylabel('Silhouette Score')\nplt.title('Silhouette Score by Cluster Count')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:18:47.325086Z","iopub.execute_input":"2022-07-14T12:18:47.326055Z","iopub.status.idle":"2022-07-14T12:20:30.828833Z","shell.execute_reply.started":"2022-07-14T12:18:47.325958Z","shell.execute_reply":"2022-07-14T12:20:30.827499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Full Models\nUsing the identified number of clustes we'll train a model a bayesian gaussian mixture model and a gaussian mixture model on the full, scaled dataset, and compare the silhouette score of the two approaches","metadata":{}},{"cell_type":"code","source":"b_gm = BayesianGaussianMixture(n_components=7, n_init=5, verbose = 0.5,tol = 0.0001, max_iter = 200).fit(scaled_df)\ngm = BayesianGaussianMixture(n_components=7, n_init=5, verbose = 0.5,tol = 0.0001, max_iter = 200).fit(scaled_df)\ngm_prediction = gm.predict(scaled_df)\nb_gm_prediction = b_gm.predict(scaled_df)\nscore_b_gm = metrics.silhouette_score(scaled_df, b_gm_prediction, metric='euclidean')\nscore_gm = metrics.silhouette_score(scaled_df, gm_prediction, metric='euclidean')\nprint(score_b_gm)\nprint(score_gm)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:22:05.586060Z","iopub.execute_input":"2022-07-14T12:22:05.586493Z","iopub.status.idle":"2022-07-14T12:27:53.154084Z","shell.execute_reply.started":"2022-07-14T12:22:05.586458Z","shell.execute_reply":"2022-07-14T12:27:53.151511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualizing Results\n\nUsing principle component analysis we can reduce our dataset to 2 dimensions and visualize the clustering below. Unfortunately we cant capture all the variance of the dataset in two dimensions as shown by plotting the explained variance of each principle component so our 2d clustering visualization isnt perfect.","metadata":{}},{"cell_type":"code","source":"gm_prediction = gm.predict(scaled_df)\npca = PCA()\npca.fit_transform(scaled_df)\npca2 = PCA(n_components=2)\npca2.fit(scaled_df)\n\n#Visualizing variance explanation of each principle component\nvariance = pca.explained_variance_\nplt.figure(figsize=(8, 6))\nplt.bar(range(len(variance)), variance, alpha=0.5, align='center', label='individual variance')\nplt.legend()\nplt.ylabel('Variance ratio')\nplt.xlabel('Principal components')\nplt.show()\n\n#Visualizing clustering\nsns.set(rc = {'figure.figsize':(8,6)})\nscaled_pca = pd.DataFrame(pca2.transform(scaled_df), columns = ['PCA1', 'PCA2'])\nscaled_pca['Group Prediction'] = gm_prediction.astype(str)\nsns.scatterplot(data = scaled_pca, x = 'PCA1', y = 'PCA2', hue = 'Group Prediction').set(title = 'Clustering Group Visualization with PCA')\n","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:27:53.641774Z","iopub.execute_input":"2022-07-14T12:27:53.642194Z","iopub.status.idle":"2022-07-14T12:28:02.102718Z","shell.execute_reply.started":"2022-07-14T12:27:53.642160Z","shell.execute_reply":"2022-07-14T12:28:02.101229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Making Prediction and Writing to File","metadata":{}},{"cell_type":"code","source":"b_gm_prediction = b_gm.predict(scaled_df)\nlabels = df['id']\nb_gm_submission = pd.DataFrame(np.array([labels, b_gm_prediction]).T,\n                                 columns = ['Id', 'Predicted'])\nb_gm_submission.to_csv('b_gm_output2.csv', index=False)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-14T12:20:30.847587Z","iopub.status.idle":"2022-07-14T12:20:30.848405Z","shell.execute_reply.started":"2022-07-14T12:20:30.847942Z","shell.execute_reply":"2022-07-14T12:20:30.847971Z"},"trusted":true},"execution_count":null,"outputs":[]}]}