{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"}],"dockerImageVersionId":30635,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"papermill":{"default_parameters":{},"duration":8.700023,"end_time":"2024-01-20T19:46:17.625939","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-01-20T19:46:08.925916","version":"2.4.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## HMS 2024 - Brain Clusters","metadata":{}},{"cell_type":"markdown","source":"Classifying different types of Hazardous Brain Activity, HBA, using EEG recordings and spectra.\nThis notebook looks at the k-means clustering of the 6-dimension probability vectors of the HBA samples.\n\nGeneral observations:
\n- The number of votes falls into two ranges: 62.6% are in 1--7 and 37.4% in 10--28.
\n- Distribution of the expert HBAs is very different between the vote ranges, most notable is that there are very few seizures in the 10--28 vote range.
\n- Is the LB test data similar to the train data? Probably not: using the train average probabilities with train data gives a KL value (e.g., CV) of 1.38 but when those probabilities are submitted as predictions the LB value is 1.07.
\n- The vote entropy gives a measure of how spread out the voting is; 89.3% have entropy at or less than a spread-over-2-bins value, and almost half have zero entropy (unanimous voting.)
\n- p-values for the hypothesis test: Ho: p_max = 0.5 and Ha: P_max > 0.5 give a measure of how convincing the votes are that the maximum votes HBA has P>0.5, i.e. is the main HBA.\n- Each patient has only about 2 HBA types used with them.
\n- Each eeg_id is from a single patient and i) mostly has one HBA type, and ii) often multiple eeg_sub_ids have the same voting values (number of votes in HBA and total votes) suggesting they were evaluated together.
\n\n\nClustering results:
\n- If clusters are determined using only the 10-28 vote range HBA samples, then 5 clusters are indicated, one for each HBA type except seizure (not surprisingly.)
\n- Using all the data, 6 clusters are indicated, one for each HBA type. Each cluster has some admixture of the others HBA types. For example the cluster center for `LRDA` is\n`[ 1.6% 6.1% 0.8% 69.9% 6.7% 15.0% ]` and shows a strong `LRDA` component, 69.9%, but also a strong 15% in the `Other` component.
\n- Adding a 7th cluster produces a new center away from the extremes: `[3.4% 12.1% 5.5% 15.7% 18.8% 44.5%]` and perhaps captures the main ambiguity in the voting.
\n- A submission could be made by using a multi-class classifier to predict the class (e.g., cluster membership) of each sample. Using 6 (or 7) clusters and assigning them perfectly to the train data gives a KL divergence score of 0.30 (0.28); this compares to a score of 0.38 if prob is set to an optimal 81.5% in the consensus HBA and 3.7% in all the others.
\n","metadata":{}},{"cell_type":"markdown","source":"## Things to use","metadata":{"papermill":{"duration":0.006481,"end_time":"2024-01-20T19:46:11.965375","exception":false,"start_time":"2024-01-20T19:46:11.958894","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n# For calculating p-values\nfrom scipy.stats import binom\n\n# For k-means\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import silhouette_score","metadata":{"papermill":{"duration":0.457861,"end_time":"2024-01-20T19:46:12.429264","exception":false,"start_time":"2024-01-20T19:46:11.971403","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:26:16.095268Z","iopub.execute_input":"2024-01-30T00:26:16.095809Z","iopub.status.idle":"2024-01-30T00:26:17.369649Z","shell.execute_reply.started":"2024-01-30T00:26:16.095733Z","shell.execute_reply":"2024-01-30T00:26:17.368262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Directory prefix for the data\nabove_dir = \"../input/hms-harmful-brain-activity-classification/\"\n# or, offline, use my local directory\n##above_dir = \"D:/Kaggle/input/hms-harmful-brain-activity-classification/\"","metadata":{"papermill":{"duration":0.013537,"end_time":"2024-01-20T19:46:12.450707","exception":false,"start_time":"2024-01-20T19:46:12.437170","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:26:17.372984Z","iopub.execute_input":"2024-01-30T00:26:17.373513Z","iopub.status.idle":"2024-01-30T00:26:17.380220Z","shell.execute_reply.started":"2024-01-30T00:26:17.373467Z","shell.execute_reply":"2024-01-30T00:26:17.378779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Functions, Etc.","metadata":{}},{"cell_type":"code","source":"# Define these since they get used a lot\nHBA_names = [\"seizure\", \"lpd\", \"gpd\", \"lrda\", \"grda\", \"other\"]\niHBA_of_expert = {\"Seizure\":0,\"LPD\":1,\"GPD\":2,\"LRDA\":3,\"GRDA\":4,\"Other\":5}\nHBA_votes = [\"seizure_vote\", \"lpd_vote\", \"gpd_vote\",\n \"lrda_vote\", \"grda_vote\", \"other_vote\"]\nHBA_probs = [\"seizure_prob\", \"lpd_prob\", \"gpd_prob\",\n \"lrda_prob\", \"grda_prob\", \"other_prob\"]\n\n# Output format for arrays\nnp.set_printoptions(precision=6, suppress=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:17.381946Z","iopub.execute_input":"2024-01-30T00:26:17.382391Z","iopub.status.idle":"2024-01-30T00:26:17.393751Z","shell.execute_reply.started":"2024-01-30T00:26:17.382351Z","shell.execute_reply":"2024-01-30T00:26:17.392082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def kld_score(solution, submission):\n '''\n Calculate the average KL divergence score.\n Ignores the \"row id\" assumed in the first column.\n '''\n sumsum = 0.0\n # Go through the probabilities\n for prob_col in solution.columns.values:\n sumsum += np.nansum(-1.0*solution[prob_col] *\n np.log(submission[prob_col] / solution[prob_col]))\n return sumsum/(len(solution))\n\n# Example of KL Divergence result from its Kaggle metric page\nsolution = pd.DataFrame({'id': range(3), 'ham': [0, 0.5, 0.5], \n 'spam': [0.1, 0.5, 0.5], 'other': [0.9, 0, 0]})\nsubmission = pd.DataFrame({'id': range(3), 'ham': [0.2, 0.3, 0.5], \n 'spam': [0.1, 0.5, 0.5], 'other': [0.7, 0.2, 0]})\n# score(solution, submission, 'id')\n# 0.160531...\n\n# Check that this simple version above gives the same value:\nkld_score(solution, submission)","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:17.395884Z","iopub.execute_input":"2024-01-30T00:26:17.396451Z","iopub.status.idle":"2024-01-30T00:26:17.424591Z","shell.execute_reply.started":"2024-01-30T00:26:17.396400Z","shell.execute_reply":"2024-01-30T00:26:17.423158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_hms_meta():\n '''\n Read in the train.csv and test.csv files.\n Add total_vote, _prob columns, and vote entropy to train_meta.\n Add extra cols to test to allow the same processing as train:\n eeg[spectro]_sub_id, eeg[spectro]_label_offset_seconds, label_id\n Make various plots of the train_meta values.\n '''\n\n # Read the test meta data\n test_meta = pd.read_csv(above_dir+\"test.csv\")\n test_meta_len = len(test_meta)\n print(\"Test has length\", test_meta_len)\n # Add columns to allow similar train/test processing:\n test_meta[\"eeg_sub_id\"] = 0\n test_meta[\"eeg_label_offset_seconds\"] = 0.0\n test_meta[\"spectrogram_sub_id\"] = 0\n test_meta[\"spectrogram_label_offset_seconds\"] = 0.0\n test_meta[\"label_id\"] = test_meta.eeg_id\n # Can decide to replace the not-real test with training data instead\n if test_meta_len > 1:\n REAL_TEST = True\n else:\n REAL_TEST = False\n print(\" --> not the real LB test data.\\n\")\n # Replace the test_meta?\n pass\n \n # Read the train meta data\n train_meta = pd.read_csv(above_dir+\"train.csv\")\n train_meta_len = len(train_meta)\n print(\"Train has length\", train_meta_len, \" with:\")\n \n # Add a total_vote column\n train_meta[\"total_vote\"] = ( train_meta[\"seizure_vote\"] +\n train_meta[\"lpd_vote\"] + train_meta[\"gpd_vote\"] +\n train_meta[\"lrda_vote\"] + train_meta[\"grda_vote\"] +\n train_meta[\"other_vote\"] )\n # Add a max_vote column (i.e. number of votes in the expert consensus)\n train_meta[\"max_vote\"] = np.max(np.array([train_meta[\"seizure_vote\"] ,\n train_meta[\"lpd_vote\"] , train_meta[\"gpd_vote\"] ,\n train_meta[\"lrda_vote\"] , train_meta[\"grda_vote\"] ,\n train_meta[\"other_vote\"]]), axis=0)\n \n # Show various unique numbers\n for this_col in [\"label_id\",\"eeg_id\",\"spectrogram_id\",\n \"patient_id\",\"total_vote\"]:\n print(\" \", len(train_meta[this_col].unique()),\n \"unique \"+this_col+\" values.\")\n \n # Look at the total votes values: 1 to 28 with missing 8 and 9\n print(\"\\nHistogram of the total votes\")\n plt.figure(figsize=(6,3))\n plt.hist(train_meta[\"total_vote\"],bins=55,log=True)\n plt.title(\"Histogram of Total Votes\")\n plt.show()\n \n # Distribution of the \"Expert consensus\" for less/more than 9 votes\n allvt = train_meta.expert_consensus.value_counts()\n less9 = train_meta[train_meta[\"total_vote\"] < 9].expert_consensus.value_counts()\n more9 = train_meta[train_meta[\"total_vote\"] > 9].expert_consensus.value_counts()\n vnot3 = train_meta[train_meta[\"total_vote\"] != 3].expert_consensus.value_counts()\n print(\"Counts for votes less than 9:\\n\",less9[allvt.index])\n print(\"\\nCounts for votes more than 9:\\n\",more9[allvt.index])\n ##print(\"\\nCounts for votes not equal to 3:\\n\",vnot3[allvt.index])\n \n # Create _prob values from the _vote values\n print(\"\\nHistograms of the probabilites of the different HBAs:\")\n print(\" (note that the large Prob=0 bin is not included.)\")\n for col_pre in HBA_names:\n train_meta[col_pre + \"_prob\"] = (train_meta[col_pre + \"_vote\"] / \n train_meta[\"total_vote\"] )\n # Show the probability histogram for each type\n plt.figure(figsize=(6,1.5))\n plt.hist(train_meta[col_pre + \"_prob\"],bins=20,range=(0.02,1))\n plt.ylim(0,len(train_meta)/5)\n plt.title(\"Histogram of \"+col_pre+\"_prob\")\n plt.show()\n \n # Calculate the entropy for each row (~ amount of vote variation)\n print(\"Calculating voting entropy values ...\")\n def calc_entropy(row):\n the_probs = np.clip(row[16:21+1].values.astype(float), 1.e-8,1.0)\n return np.nansum(the_probs * -1*np.log(the_probs))\n # Add an entropy column\n train_meta[\"entropy\"] = train_meta.apply(calc_entropy, axis=1)\n \n return train_meta, test_meta\n","metadata":{"papermill":{"duration":0.021479,"end_time":"2024-01-20T19:46:12.478233","exception":false,"start_time":"2024-01-20T19:46:12.456754","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:26:17.428027Z","iopub.execute_input":"2024-01-30T00:26:17.428640Z","iopub.status.idle":"2024-01-30T00:26:17.449989Z","shell.execute_reply.started":"2024-01-30T00:26:17.428603Z","shell.execute_reply":"2024-01-30T00:26:17.448749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prob_prob_scatter(name1, name2, probs2plot, clust_ids):\n '''\n Make a prob1 vs prob2 scatter plot.\n Include an x at the cluster centers in chosen axes.\n Use sqrt scaling.\n External: HBA_probs, iHBA_of_expert[ ], clust_probs\n '''\n # Color-code the vectors by their k-means label, arbitrary\n kmclrs = 2*[\"red\",\"blue\",\"green\",\"black\",\"purple\",\"orange\"]\n clstclrs = []\n for ilab in clust_ids:\n clstclrs.append(kmclrs[ilab])\n \n ixax = iHBA_of_expert[name1]\n iyax = iHBA_of_expert[name2]\n lenprob = len(probs2plot)\n plt.figure(figsize=(5,5))\n plt.scatter(np.sqrt(probs2plot[HBA_probs[ixax]]) + \n 0.04*(np.random.rand(lenprob)-0.5),\n np.sqrt(probs2plot[HBA_probs[iyax]]) + \n 0.04*(np.random.rand(lenprob)-0.5),\n s=3, c=clstclrs, alpha=0.02)\n # Add the centers\n for iclust in range(0,len(clust_probs)):\n plt.plot(np.sqrt([clust_probs[iclust,ixax]]),\n np.sqrt([clust_probs[iclust,iyax]]),\n c=kmclrs[iclust],marker=\"x\",markersize=15)\n plt.xlabel(\"sqrt( \"+name1+\" )\")\n plt.ylabel(\"sqrt( \"+name2+\" )\")\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:17.451874Z","iopub.execute_input":"2024-01-30T00:26:17.452282Z","iopub.status.idle":"2024-01-30T00:26:17.468725Z","shell.execute_reply.started":"2024-01-30T00:26:17.452247Z","shell.execute_reply":"2024-01-30T00:26:17.467353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"## Get and look at the csv meta data","metadata":{}},{"cell_type":"code","source":"# Read in the meta data, routine also looks at the train values\ntrain_meta, test_meta = read_hms_meta()","metadata":{"papermill":{"duration":2.42584,"end_time":"2024-01-20T19:46:14.909979","exception":false,"start_time":"2024-01-20T19:46:12.484139","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:26:17.470224Z","iopub.execute_input":"2024-01-30T00:26:17.470584Z","iopub.status.idle":"2024-01-30T00:26:33.833192Z","shell.execute_reply.started":"2024-01-30T00:26:17.470551Z","shell.execute_reply":"2024-01-30T00:26:33.831648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The columns in train_meta\ntrain_meta.info()","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:33.835003Z","iopub.execute_input":"2024-01-30T00:26:33.836090Z","iopub.status.idle":"2024-01-30T00:26:33.867939Z","shell.execute_reply.started":"2024-01-30T00:26:33.836050Z","shell.execute_reply":"2024-01-30T00:26:33.866761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The columns in test_meta\ntest_meta.info()","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:33.869533Z","iopub.execute_input":"2024-01-30T00:26:33.870110Z","iopub.status.idle":"2024-01-30T00:26:33.885961Z","shell.execute_reply.started":"2024-01-30T00:26:33.870074Z","shell.execute_reply":"2024-01-30T00:26:33.884271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Look at the vote entropies","metadata":{}},{"cell_type":"code","source":"if True:\n print(\"Histogram of the votes entropies\")\n plt.figure(figsize=(6,3))\n plt.hist(train_meta[\"entropy\"],bins=50,log=True)\n plt.title(\"Histogram of Vote Entropy (~ vote variation)\")\n plt.show()\n \n # Show the Vote Entropy vs the Number of Votes\n plt.figure(figsize=(8,6))\n plt.scatter(train_meta['total_vote'] + 0.7*(np.random.rand(len(train_meta))-0.5),\n train_meta['entropy'] + 0.05*(np.random.rand(len(train_meta))-0.5),\n s=3, alpha=0.02)\n # For reference, plot entropy values if spread evenly into n bins\n ref_ents = []\n for ispread in [1,2,3,4,5,6]:\n spread_ent = np.log(ispread)\n ref_ents.append(spread_ent)\n plt.plot([ispread,28],[spread_ent,spread_ent],\n lw=2, c='pink', alpha=0.5)\n plt.text(24.0, spread_ent+0.03, \"{} x p=1/{}\".format(\n ispread,ispread))\n plt.title(\"Vote Entropy vs Number of Votes\") \n plt.ylim(-0.05,1.90)\n plt.ylabel(\"Entropy of the Votes\")\n plt.xlabel(\"Number of Votes\")\n plt.show()\n # List the reference entropy values\n fmtstr = (len(ref_ents)-1) * '{:.4f}, ' + '{:.4f}'\n print(\" The entropy reference lines are at:\",\n fmtstr.format(*ref_ents),\"\\n\")\n frac_z = sum(train_meta[\"entropy\"] < 0.02)/len(train_meta[\"entropy\"])\n frac_lt2 = sum(train_meta[\"entropy\"] < 0.70)/len(train_meta[\"entropy\"])\n frac_lt3 = sum(train_meta[\"entropy\"] < 1.10)/len(train_meta[\"entropy\"])\n print(\" The fraction with zero entropy is {:.1f}%\".format(\n 100.0*frac_z))\n print(\" The fraction below 2 x p=1/2 is {:.1f}%\".format(\n 100.0*frac_lt2))\n print(\" The fraction below 3 x p=1/3 is {:.1f}%\".format(\n 100.0*frac_lt3))\n ","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:33.887886Z","iopub.execute_input":"2024-01-30T00:26:33.888274Z","iopub.status.idle":"2024-01-30T00:26:35.221616Z","shell.execute_reply.started":"2024-01-30T00:26:33.888240Z","shell.execute_reply":"2024-01-30T00:26:35.220469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### p-values for votes","metadata":{}},{"cell_type":"code","source":"# p-value for rejecting Ho: p_max = 0.5\n# and so accepting Ha: p_max > 0.5\n## x_votes = 12\n## n_votes = 15\n## 1 - binom.cdf(x_votes-1, n_votes, 0.5)\n\n# Get p-values for all the rows\np_vals = []\nfor (x_votes, n_votes) in zip(train_meta.max_vote, train_meta.total_vote):\n if x_votes == n_votes:\n # If it is unanimous assume very high confidence, very low p-value\n p_vals.append(np.exp(-7.2))\n else:\n # Assign the p-value based on binomial distribution \n p_vals.append(np.clip((1 - binom.cdf(x_votes-1, n_votes, 0.5)),0.001,1.0))\n \np_vals = np.array(p_vals)\nlog_p_vals = np.log(p_vals)","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:35.223913Z","iopub.execute_input":"2024-01-30T00:26:35.224369Z","iopub.status.idle":"2024-01-30T00:26:47.265922Z","shell.execute_reply.started":"2024-01-30T00:26:35.224327Z","shell.execute_reply":"2024-01-30T00:26:47.264557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Various plots using the p-values\np_markers = [0.001,0.01,0.05,0.20]\n\nplt.figure(figsize=(6,5))\nplt.scatter(log_p_vals + 0.12*(np.random.rand(len(train_meta))-0.5),\n train_meta.entropy + 0.03*(np.random.rand(len(train_meta))-0.5),\n s=3,alpha=0.02)\nfor alpha in p_markers:\n plt.plot(2*[np.log(alpha)],[0.0,1.4],c='orange',lw=1)\n plt.text(np.log(alpha)-0.2, 1.43,str(100*alpha)+\"%\",c='orange')\nplt.xlabel(\"log( p-value )\")\nplt.ylabel(\"One HBA < - - Entropy - - > Spread among HBAs\")\nplt.title(\"Vote Entropy vs p-value\")\nplt.show()\n\nplt.figure(figsize=(6,5))\nplt.scatter(log_p_vals + 0.18*(np.random.rand(len(train_meta))-0.5),\n train_meta.max_vote/train_meta.total_vote +\n 0.03*(np.random.rand(len(train_meta))-0.5),\n s=3,alpha=0.02)\nfor alpha in p_markers:\n plt.plot(2*[np.log(alpha)],[0.01,0.95],c='orange',lw=1)\n plt.text(np.log(alpha)-0.2, 0.98,str(100*alpha)+\"%\",c='orange')\nplt.xlabel(\"log( p-value )\")\nplt.ylabel(\"Votes in max HBA / Total votes\")\nplt.ylim(0.2,1.05) # 0.2 is minimum possible\nplt.title(\"Max Votes / Total Votes vs p-value\")\nplt.show()\n\n\nprint(\"\\nFraction with p-value < 1.0% is {:.1f}%\\n\".format(\n 100.0 * sum(p_vals < 0.01) / len(p_vals)))\nprint(\"Fraction with p-value > 20.0% is {:.1f}%\\n\".format(\n 100.0 * sum(p_vals > 0.20) / len(p_vals)))","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:47.267571Z","iopub.execute_input":"2024-01-30T00:26:47.268138Z","iopub.status.idle":"2024-01-30T00:26:48.250349Z","shell.execute_reply.started":"2024-01-30T00:26:47.268082Z","shell.execute_reply":"2024-01-30T00:26:48.248714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Identical voting within an eeg\nBy using the number of rows when grouped in different ways, we can look into @patrob 's discussion comment about [multiple identical voting distributions](https://www.kaggle.com/competitions/hms-harmful-brain-activity-classification/discussion/470645#2618899).\n","metadata":{}},{"cell_type":"code","source":"# Look into the idea that a given patient_id--eeg_id will have sub ids\n# with mostly similar expert consensus and vote numbers.\n\nif True: # results are in comments below\n \n # To show what happens if there is no similarity,\n # include randomly assigned HBA to all rows -- using 2 or 3 choices.\n train_meta[\"rand_id\"] = np.random.choice(\n 3, size=len(train_meta), replace=True, p=None)\n\n grouped = train_meta[[\"eeg_id\",\"eeg_sub_id\",\"patient_id\",\"spectrogram_id\",\n \"expert_consensus\",\"rand_id\",\n \"total_vote\",\"max_vote\"]].groupby(\n [\"patient_id\",\"eeg_id\",\"expert_consensus\",\n \"total_vote\",\"max_vote\"]).count()\n # Histogram of \"number of identicals\" in eegs\n plt.figure(figsize=(5,3))\n plt.hist(grouped[grouped.eeg_sub_id < 50].rand_id,bins=49,log=True)\n plt.xlabel(\"Number of identical votes in the eeg\")\n plt.ylabel(\"Number of eegs\")\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:48.252078Z","iopub.execute_input":"2024-01-30T00:26:48.252508Z","iopub.status.idle":"2024-01-30T00:26:48.928837Z","shell.execute_reply.started":"2024-01-30T00:26:48.252470Z","shell.execute_reply":"2024-01-30T00:26:48.927472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Number of types of HBAs for a patient: looks like about 2 for each.\n# 1950 rows, grouped by: [\"patient_id\"]\n# 3625 rows, grouped by: [\"patient_id\",\"expert_consensus\"]\n# 3799 rows, grouped by: [\"patient_id\",\"rand_id\"] <-- 2 coices, 0,1\n# 5533 rows, grouped by: [\"patient_id\",\"rand_id\"] <-- 3 choices, 0,1,2\n\n# Variety of HBAs and votes in patient--eeg combinations\n# 17089 rows, grouped by: [\"patient_id\",\"eeg_id\"]\n# 18013 rows, grouped by: [\"patient_id\",\"eeg_id\",\"expert_consensus\"]\n# 19783 rows, grouped by: [\"patient_id\", . . . _consensus\",\"total_vote\"]\n# 20072 rows, grouped by: [\"patient_id\", . . . _consensus\",\"total_vote\",\"max_vote\"]\n# 26266 rows, grouped by: [\"patient_id\",\"eeg_id\",\"rand_id\"] <-- 2 coices, 0,1\n \n\n# Show the dataframe. \n# Can select by the number of \"identicals\" = the count is in non-groupby columns\ngrouped[grouped.eeg_sub_id > 5]","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:48.934037Z","iopub.execute_input":"2024-01-30T00:26:48.934547Z","iopub.status.idle":"2024-01-30T00:26:48.974985Z","shell.execute_reply.started":"2024-01-30T00:26:48.934492Z","shell.execute_reply":"2024-01-30T00:26:48.973593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"## Make a constant-probabilities Submission\nMake a submission file with constant probabilities: all 1/6, the mean probabilities, etc. Evaluate those probabilites vs the train probabilites using the KL Divergence metric. Can compare those KL values with the LB values to see how different train and test may be.","metadata":{"papermill":{"duration":0.009196,"end_time":"2024-01-20T19:46:15.183204","exception":false,"start_time":"2024-01-20T19:46:15.174008","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Create a 'solution' of actual probabilities from the training data\nsolution_train = train_meta[[\"eeg_id\"] + HBA_votes]\n# In solution_train replace the _votes values with probabilities\nfor col_pre in HBA_names:\n solution_train.loc[:, col_pre + \"_vote\"] = train_meta[col_pre + \"_prob\"]\n \n##solution_train","metadata":{"papermill":{"duration":0.031041,"end_time":"2024-01-20T19:46:15.223836","exception":false,"start_time":"2024-01-20T19:46:15.192795","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:26:48.976903Z","iopub.execute_input":"2024-01-30T00:26:48.977742Z","iopub.status.idle":"2024-01-30T00:26:49.000930Z","shell.execute_reply.started":"2024-01-30T00:26:48.977693Z","shell.execute_reply":"2024-01-30T00:26:48.999128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make a constant-probabilites 'submission' dataframe\nsubmission_train = solution_train.copy()\n\n# Calculate the mean prob values\nmean_probs = solution_train.iloc[ : , 1:].mean().values\nprint(\"Mean train probabilities:\", mean_probs)\n\n# Put values into the 'submission' dataframe\n##submit_probs = 6*[1/6]\nsubmit_probs = mean_probs\n# excentuate/reduce the differences between the mean probs\n##expon = 0.25\n##submit_probs = mean_probs**expon / np.sum(mean_probs**expon)\n#\nprint(\"Submitting these prob.s:\", submit_probs)\nfor iprob, col_pre in enumerate(HBA_names):\n submission_train.loc[ : , col_pre + \"_vote\"] = submit_probs[iprob]\n\n##submission_train","metadata":{"papermill":{"duration":0.037395,"end_time":"2024-01-20T19:46:15.271065","exception":false,"start_time":"2024-01-20T19:46:15.233670","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:29:14.579215Z","iopub.execute_input":"2024-01-30T00:29:14.579784Z","iopub.status.idle":"2024-01-30T00:29:14.608020Z","shell.execute_reply.started":"2024-01-30T00:29:14.579728Z","shell.execute_reply":"2024-01-30T00:29:14.606619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Run the KL divergence metric\nkld_score(solution_train, submission_train)\n\n# A small improvement using the mean probs:\n# train KL LB KL\n# A constant 1/6 : 1.4023 1.09 v11\n# Mean probs^.25 : 1.3926 1.08 v10\n# Mean probs^.50 : 1.3856 1.07 v9\n# Mean probs^.75 : 1.3815\n# The mean probs : * 1.3801 1.07 v3\n# Mean probs^1.25 : 1.3815\n# Mean probs^1.5 : 1.3855 1.07 v6\n# Mean probs^2.0 : 1.4016 1.09 v8\n# * = minimum, very broad","metadata":{"papermill":{"duration":0.017153,"end_time":"2024-01-20T19:46:15.297674","exception":false,"start_time":"2024-01-20T19:46:15.280521","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:29:19.720229Z","iopub.execute_input":"2024-01-30T00:29:19.720705Z","iopub.status.idle":"2024-01-30T00:29:19.748940Z","shell.execute_reply.started":"2024-01-30T00:29:19.720666Z","shell.execute_reply":"2024-01-30T00:29:19.747249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Assemble the submission using the desired prob.s\n# Start with a dataframe with just the eeg_id column from test_meta\ntest_submit = test_meta[[\"eeg_id\"]].copy()\n# Add \"_vote\" columns with the predicted probabilities\nfor iprob, col_pre in enumerate(HBA_names):\n test_submit[col_pre + \"_vote\"] = submit_probs[iprob]\n\nprint(test_submit)\n\n# Output the file\ntest_submit.to_csv(\"submission.csv\", header=True, \n index=False, na_rep='', float_format='%.6f')\n# Look at the file\n##!more submission.csv\n","metadata":{"papermill":{"duration":0.029538,"end_time":"2024-01-20T19:46:15.354930","exception":false,"start_time":"2024-01-20T19:46:15.325392","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-30T00:29:28.163644Z","iopub.execute_input":"2024-01-30T00:29:28.164182Z","iopub.status.idle":"2024-01-30T00:29:28.185185Z","shell.execute_reply.started":"2024-01-30T00:29:28.164134Z","shell.execute_reply":"2024-01-30T00:29:28.183580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
","metadata":{"papermill":{"duration":0.009263,"end_time":"2024-01-20T19:46:15.433291","exception":false,"start_time":"2024-01-20T19:46:15.424028","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"## Clustering of HBA Probability Vectors","metadata":{}},{"cell_type":"code","source":"# Fractions of HBAs in regions of the votes--entropy plane\n#\n# votes<9 is 62.6% , just votes=3 is 48.6%\n# votes<9 and entropy<0.01 is 44.5%\n# votes<9 and entropy<0.70 is 59.5%\n#\n# votes>9 is 37.4% , just votes=15 is 10.0%\n# votes>9 and entropy<0.01 is 3.3% <-- very few are unanimous\n# votes>9 and entropy<0.70 is 19.4%\n\n##print(len(train_meta[(train_meta.total_vote > 9) &\n## (train_meta.entropy < 0.70)]) / len(train_meta))\n\n# Select all or a subset of the HBA samples to make clusters...\n# Down-select to votes>9 ?\n##prob_vectors = train_meta.loc[(train_meta.total_vote > 9), HBA_probs]\n# or... exclude votes=3 ?\n##prob_vectors = train_meta.loc[(train_meta.total_vote != 3), HBA_probs]\n# or... Use all HBAs ?\nprob_vectors = train_meta.loc[(train_meta.total_vote > -1), HBA_probs]\n\n##prob_vectors","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:49.090650Z","iopub.execute_input":"2024-01-30T00:26:49.091208Z","iopub.status.idle":"2024-01-30T00:26:49.132113Z","shell.execute_reply.started":"2024-01-30T00:26:49.091158Z","shell.execute_reply":"2024-01-30T00:26:49.130584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Mean and standard deviation of each HBA's prob.s\nprint(prob_vectors.apply(np.mean,axis=0).values)\nprint(prob_vectors.apply(np.std,axis=0).values)\n\n# For >9 votes: The Seizures are under-represented\n# [0.05205155 0.17736811 0.20306643 0.13129067 0.12840864 0.30781461]\n# [0.11620024 0.28729203 0.31122879 0.22205968 0.22739241 0.30195994]\n# Excluding votes=3 Some more Seizures\n# [0.07298742 0.19346076 0.19402469 0.11244978 0.10907672 0.31800064]\n# [0.18098827 0.31381592 0.31565785 0.21338437 0.21958587 0.32906653]\n# For all HBAs:\n# [0.20831852 0.13211966 0.12853332 0.13891264 0.1792938 0.21282206]\n# [0.37827328 0.27772958 0.27617022 0.28005733 0.33636886 0.31519539]","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:49.134274Z","iopub.execute_input":"2024-01-30T00:26:49.134848Z","iopub.status.idle":"2024-01-30T00:26:49.159220Z","shell.execute_reply.started":"2024-01-30T00:26:49.134785Z","shell.execute_reply":"2024-01-30T00:26:49.157679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use the k-means routine in sklearn:\n# KMeans(n_clusters=8, *, init='k-means++', n_init='auto', max_iter=300,\n# tol=0.0001, verbose=0, random_state=None, copy_x=True, algorithm='lloyd')\n\n# Determine the appropriate number of clusters\nif True:\n # Look for the \"elbow\" number of clusters\n # Use a fraction of the samples\n prob_array = np.array(prob_vectors[0::10])\n maxclsts = 15\n inertias = []\n silhos = []\n for iclust in range(2, maxclsts+1):\n kmeans = KMeans(n_clusters = iclust, init = 'k-means++', n_init = 10,\n max_iter = 300, random_state=None)\n kmeans.fit(prob_array)\n inertias.append(kmeans.inertia_)\n # This metric has \"The best value is 1 and the worst value is -1\"\n silhos.append(silhouette_score(prob_array, kmeans.labels_))\n # Plot the metrics vs number of clusters\n plt.plot(range(2, maxclsts+1), inertias, 'b-')\n plt.plot(range(2, maxclsts+1), max(inertias)*np.array(silhos), 'g-')\n plt.title('Clustering Metrics vs Number of Clusters')\n plt.xlabel('Number of clusters')\n plt.ylabel('inertia (blue), scaled silhouette (green)')\n plt.show()\n\n# For Votes>9: The inertia elbow (blue) is at 5, peak of silho' is also at 5.\n# Not=3: The inertia elbow (blue) is at 6, peak of silho' is also at 6.\n# Using all: The inertia elbow (blue) is at 6, peak of silho' is also at 6.","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:26:49.161949Z","iopub.execute_input":"2024-01-30T00:26:49.162561Z","iopub.status.idle":"2024-01-30T00:27:31.744832Z","shell.execute_reply.started":"2024-01-30T00:26:49.162491Z","shell.execute_reply":"2024-01-30T00:27:31.743439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find the centers for the optimum (or other) number of clusters \n##iclust = 5 # when using votes > 9\niclust = 6+1 # when using all the vectors - adding a 7th cluster\n\nprob_array = np.array(prob_vectors)\nkmeans = KMeans(n_clusters = iclust, init = 'k-means++', n_init = 10,\n max_iter = 300, random_state=None)\nkmeans.fit(prob_array)\n# Can get the cluster center coord.s/probs from:\nclust_probs = kmeans.cluster_centers_\nprint(\"The centers for {} clusters:\".format(iclust))\nprint(clust_probs)\n\n# Using votes > 9 rows only\n# The centers for 5 clusters: * ordered by main column *\n# [ none that is mostly seizure ]\n# [0.03971239 0.75258142 0.03044835 0.05662892 0.00803234 0.11259658]\n# [0.12050526 0.04732145 0.70046799 0.00519864 0.02985535 0.09665131]\n# [0.04379057 0.11936496 0.01280059 0.53104421 0.05826677 0.2347329 ]\n# [0.00887504 0.02846734 0.05263628 0.04578226 0.60112928 0.2631098 ]\n# [0.02070779 0.05335455 0.03500729 0.04696622 0.06429828 0.77966588]\n\n# Using all except votes=3\n# The centers for 6 clusters: * ordered by main column *\n# [0.73180621 0.0631708 0.04338586 0.02012522 0.01133559 0.13017633]\n# [0.0298963 0.79677792 0.02293553 0.04389117 0.00579372 0.10070537]\n# [0.08950668 0.05180871 0.71981674 0.00475706 0.02952677 0.10458404]\n# [0.03100587 0.11865955 0.01128878 0.53514721 0.05215337 0.25174522]\n# [0.00691587 0.02549976 0.05270252 0.04168656 0.61507058 0.25812471]\n# [0.01257443 0.04000609 0.02475847 0.03336829 0.05108967 0.83820305]\n\n# Using All rows\n# The centers for 6 clusters: * ordered by main column * \n# [0.97144241 0.00613458 0.00485837 0.00248338 0.00133784 0.01374342]\n# [0.03762207 0.78874927 0.02092333 0.04629211 0.00839672 0.0980165 ]\n# [0.08720708 0.05977927 0.72124703 0.00427439 0.0329343 0.09455794]\n# [0.01610773 0.06060393 0.00790338 0.69877109 0.06680986 0.149804 ]\n# [0.0033314 0.00741154 0.01295589 0.02588149 0.88386741 0.06655228]\n# [0.01802665 0.03991739 0.02326481 0.04377798 0.07160484 0.80340834]\n#\n# The centers for 7 clusters: * ordered by main column *\n# [0.97295318 0.00590648 0.00483242 0.00223349 0.00114005 0.01293438]\n# [0.0359414 0.8263946 0.02026326 0.0416359 0.0039719 0.07179293]\n# [0.09039322 0.05736363 0.73959869 0.0036609 0.02577728 0.08320629]\n# [0.01311906 0.05057075 0.00535362 0.76954905 0.05295756 0.10844996]\n# [0.00294371 0.00356274 0.00838875 0.02236721 0.92515794 0.03757965]\n# [0.00794713 0.01162694 0.01349516 0.01482245 0.02306207 0.92904625]\n# The new 7th cluster:\n# [0.03401079 0.12059055 0.05457865 0.15692139 0.18847755 0.44542107] ","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:31.747055Z","iopub.execute_input":"2024-01-30T00:27:31.747664Z","iopub.status.idle":"2024-01-30T00:27:33.481051Z","shell.execute_reply.started":"2024-01-30T00:27:31.747617Z","shell.execute_reply":"2024-01-30T00:27:33.479560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Not sure what is the best way to view these...\n# Try pair scatter plots with sqrt scaling\n# Color code by cluster id\nclust_ids = kmeans.labels_\n\nprob_prob_scatter(\"LPD\",\"GRDA\", prob_vectors, clust_ids)\n\nprob_prob_scatter(\"LRDA\",\"Other\", prob_vectors, clust_ids)\n\nprob_prob_scatter(\"Seizure\",\"GPD\", prob_vectors, clust_ids)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:33.482807Z","iopub.execute_input":"2024-01-30T00:27:33.483276Z","iopub.status.idle":"2024-01-30T00:27:40.175968Z","shell.execute_reply.started":"2024-01-30T00:27:33.483241Z","shell.execute_reply":"2024-01-30T00:27:40.174492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluate KL metric when cluster probabilites are assigned","metadata":{}},{"cell_type":"code","source":"# Assign clusters to the train_meta\ntrain_meta[\"clust_id\"] = kmeans.predict(np.array(train_meta[HBA_probs]))\nclust_probs = kmeans.cluster_centers_\n\n\n# For comparison:\n# Use the expert consensus classes with optimal near-'unit' probabilites:\nif False:\n for irow in train_meta.index:\n train_meta.loc[irow,\"clust_id\"] = iHBA_of_expert[\n train_meta.loc[irow,\"expert_consensus\"]]\n pyes = 0.815\n pnot = (1-pyes)/5.0\n clust_probs = np.array([[pyes-pnot,0.0,0.0,0.0,0.0,0.0],\n [0.0,pyes-pnot,0.0,0.0,0.0,0.0],\n [0.0,0.0,pyes-pnot,0.0,0.0,0.0],\n [0.0,0.0,0.0,pyes-pnot,0.0,0.0],\n [0.0,0.0,0.0,0.0,pyes-pnot,0.0],\n [0.0,0.0,0.0,0.0,0.0,pyes-pnot]]) + pnot\n\n\n# The cluster probability vectors\nprint(clust_probs)\n\n# main component of each cluster\nmax_probs = np.argmax(clust_probs,axis=1)\nclust_names = []\nfor iclust in range(len(clust_probs)):\n clust_names.append(HBA_names[max_probs[iclust]])\nclust_names","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:40.177929Z","iopub.execute_input":"2024-01-30T00:27:40.178474Z","iopub.status.idle":"2024-01-30T00:27:40.223304Z","shell.execute_reply.started":"2024-01-30T00:27:40.178426Z","shell.execute_reply":"2024-01-30T00:27:40.222088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compare cluster counts and expert counts\nprint(train_meta[\"clust_id\"].value_counts())\nprint(train_meta[\"expert_consensus\"].value_counts())\n\n# These agree that the most are in Seizure and the Least are in LPD.","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:40.227063Z","iopub.execute_input":"2024-01-30T00:27:40.227796Z","iopub.status.idle":"2024-01-30T00:27:40.255295Z","shell.execute_reply.started":"2024-01-30T00:27:40.227738Z","shell.execute_reply":"2024-01-30T00:27:40.253499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Start with the previous submission_train and replace the const probs\n# with the appropriate cluster probs: \n# Go through the 6 probability columns and\n# use cluster id to select the correct cluster prob value\nfor iprob in range(len(clust_probs[0])):\n this_col_probs = clust_probs[: , iprob]\n submission_train[HBA_votes[iprob]] = this_col_probs[train_meta[\"clust_id\"]]\n\nsubmission_train","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:40.257324Z","iopub.execute_input":"2024-01-30T00:27:40.258266Z","iopub.status.idle":"2024-01-30T00:27:40.290977Z","shell.execute_reply.started":"2024-01-30T00:27:40.258208Z","shell.execute_reply":"2024-01-30T00:27:40.289465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Evaluate the KL divergence when predictions are the cluster probs\nkld_score(solution_train, submission_train)\n\n# Using the \"votes>9\" 5 clusters : 0.8852\n# Using expert consensus, 90%,2% : 0.4163\n# \" 81.5% 3.7% : 0.3841\n# Using votes-not-3 6 clusters : 0.3705\n# Using the \"all rows\" 6 clusters : 0.3001\n# Using all rows and 6+1 clusters : 0.2780 (+/-0.0020, variation in 7th cluster)\n\n# This suggests that a multi-class classifier using the clusters could do well. ","metadata":{"execution":{"iopub.status.busy":"2024-01-30T00:27:40.292993Z","iopub.execute_input":"2024-01-30T00:27:40.294139Z","iopub.status.idle":"2024-01-30T00:27:40.326412Z","shell.execute_reply.started":"2024-01-30T00:27:40.294078Z","shell.execute_reply":"2024-01-30T00:27:40.325082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"
","metadata":{}}]}