{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# let's set things up\r\n",
    "from IPython.core.interactiveshell import InteractiveShell\r\n",
    "import matplotlib as mpl\r\n",
    "import matplotlib.pyplot as plt\r\n",
    "import numpy as np\r\n",
    "import pandas as pd\r\n",
    "import seaborn as sns\r\n",
    "from sklearn.datasets import load_diabetes\r\n",
    "\r\n",
    "InteractiveShell.ast_node_interactivity = \"all\"\r\n",
    "%matplotlib inline\r\n",
    "plt.style.use('default')\r\n",
    "sns.set()\r\n",
    "pd.options.display.float_format = '{:,.2f}'.format"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# original data source: https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt\r\n",
    "db_loc = \"./data/diabetes.tab.txt\"\r\n",
    "db_data = pd.read_csv(db_loc, sep='\\t', header=(0))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "print(type(db_data))\r\n",
    "db_data.head()\r\n",
    "db_data.describe()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "db_data.dtypes\r\n",
    "db_data.info()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# let's have a look at the target data\r\n",
    "g = sns.displot(data=db_data, x='Y', kde=True, rug=True)\r\n",
    "#g.set_title('Disease Progression Distribution');"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# and lets look at that by sex, not that I know which is which\r\n",
    "sns.displot(data=db_data, x=\"Y\", kde=True, col=\"SEX\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Similar shapes for both sexes. But there seems to be a higher percentage high values for sex = 2.\r\n",
    "\r\n",
    "I'd like to look at the influence of BMI on the progression value. But, the above approach won't work given the continuous, rather than categorical, nature of the BMI feature. So, let's look at BMI by sex first."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "sns.displot(data=db_data, x=\"BMI\", kde=True, col=\"SEX\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "So, what do you think, are males 1 or 2? I'm guessing 2, but...\r\n",
    "\r\n",
    "Will need to see if I can find any reliable information on the dataset. Seems to be scarce."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# how about sex and blood pressure\r\n",
    "sns.displot(data=db_data, x=\"BP\", kde=True, col=\"SEX\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Well that almost looks like the reverse of the previous chart.\r\n",
    "\r\n",
    "Okay let's perhaps look at the distributions of all the features."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "fig=plt.figure(figsize=(20,20))\r\n",
    "for i,col in enumerate(db_data.drop(['Y'],axis=1)):\r\n",
    "    ax=fig.add_subplot(5,2,i+1);\r\n",
    "    sns.histplot(db_data[col]);"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Most of the distributions appear to be somewhat skewed. Which is probably why the scikit-learn data was mean centered and scaled by the standard deviation times the number of samples (i.e. the sum of squares of each column totals 1).\r\n",
    "\r\n",
    "There appears to only be one categorical feature, SEX. Though the distribution for S4 (tch) looks a bit odd to me.\r\n",
    "\r\n",
    "Let's have a look at the box plots for the continuous variables — after a bit of a side step."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# count by sex\r\n",
    "c_sx = db_data.SEX.count()\r\n",
    "c_s1 = db_data[db_data['SEX'] == 1.0].SEX.count()\r\n",
    "print(f\"sex 1: {c_s1}, sex 2: {c_sx - c_s1}\\n\")\r\n",
    "\r\n",
    "s4 = db_data.S4.unique()\r\n",
    "# sorted(s4)\r\n",
    "print(f's4 unique: {len(s4)}')\r\n",
    "# The American Board of Internal Medicine uses an adult TSH reference range of 0.5–4.0 mIU/L.\r\n",
    "# how many values are outside this range\r\n",
    "s4_low = db_data[db_data['S4'] < 0.5]\r\n",
    "s4_high = db_data[db_data['S4'] > 4.0]\r\n",
    "c_s4_h_all = len(s4_high)\r\n",
    "print(f's4_low: {len(s4_low)}, s4 ok: {c_sx - len(s4_low) - len(s4_high)}, s4_high: {c_s4_h_all}\\n')\r\n",
    "\r\n",
    "# high by sex\r\n",
    "#s4_h_1 = db_data[db_data['S4'] > 4.0 and db_data['SEX'] == 1.0]\r\n",
    "s4_h_1 = db_data.query(\"S4 > 4.0 and SEX == 1.0\")\r\n",
    "c_s4_h_1 = len(s4_h_1)\r\n",
    "print(f\"s4 high sex 1: {c_s4_h_1}; s4 high sex 2: {c_s4_h_all - c_s4_h_1}\")\r\n",
    "\r\n",
    "# maybe look at turning this into a categorical value? low, ok, high or -1, 0, 1"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "boxp = pd.DataFrame(data=db_data, columns=['AGE', 'BMI', 'S3'])\r\n",
    "sns.boxplot(x=\"variable\", y=\"value\", data=pd.melt(boxp), showmeans=True);"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "boxp2 = pd.DataFrame(data=db_data, columns=['BP', 'S1', 'S2', 'S6'])\r\n",
    "sns.boxplot(x=\"value\", y=\"variable\", data=pd.melt(boxp2), showmeans=True, orient=\"h\", palette=\"Set3\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "boxp3 = pd.DataFrame(data=db_data, columns=['S4', 'S5'])\r\n",
    "sns.boxplot(x=\"value\", y=\"variable\", data=pd.melt(boxp3), showmeans=True, palette=\"Set2\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "I kind of like the horizontal display. But, either view does the job. I originally plotted all the variables on one chart. But the features with small values and/or ranges were rather compressed. So I split them up.\r\n",
    "\r\n",
    "Clearly some outliers in the data for a few of the features. But on the other hand, there doesn't appear to be a significant amount of skewness in most of the feature distributions.\r\n",
    "\r\n",
    "Now, let's look at a similar visual: the violin plot."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# display boxplot and violin plot side by side for comparison\r\n",
    "#fig, axs = plt.subplots(3, 2, figsize=(16,16), sharey='row')\r\n",
    "fig, axs = plt.subplots(3, 2, figsize=(16,16))\r\n",
    "#axs[0,1].ylim(0, 120)\r\n",
    "sns.violinplot(x=\"variable\", y=\"value\", data=pd.melt(boxp), showmeans=True, ax=axs[0,0]);\r\n",
    "# yticks = axs[0, 0].get_yticks()\r\n",
    "# print(yticks)\r\n",
    "# axs[0, 1].set_yticks(yticks)\r\n",
    "sns.boxplot(x=\"variable\", y=\"value\", data=pd.melt(boxp), showmeans=True, ax=axs[0,1]);\r\n",
    "axs[0, 1].sharey(axs[0, 0])\r\n",
    "# spent a bunch of time trying to get y-axis on top row to have same ylim and to show tick values\r\n",
    "#   turned out to be simpler than I thought, also made x-axis on lower rows match up\r\n",
    "#   don't know if this could be done someway easier in Seaborn, don't think so\r\n",
    "#plt.setp(axs[0,1], xticklabels=[])\r\n",
    "# axs[0, 0].tick_params(axis='both', which='both', labelsize=7, labelbottom=True)\r\n",
    "# axs[0, 1].tick_params(axis='both', which='both', left=True, labelsize=7)\r\n",
    "# for tick in axs[0, 1].get_yticklabels():\r\n",
    "#     tick.set_visible(True)\r\n",
    "\r\n",
    "sns.violinplot(x=\"value\", y=\"variable\", data=pd.melt(boxp2), showmeans=True, orient=\"h\", palette=\"Set3\", ax=axs[1,0]);\r\n",
    "l1, r1 = axs[1, 0].get_xlim()\r\n",
    "axs[1, 1].set_xlim(left=l1, right=r1)\r\n",
    "sns.boxplot(x=\"value\", y=\"variable\", data=pd.melt(boxp2), showmeans=True, orient=\"h\", palette=\"Set3\", ax=axs[1,1]);\r\n",
    "\r\n",
    "sns.violinplot(x=\"value\", y=\"variable\", data=pd.melt(boxp3), showmeans=True, palette=\"Set2\", ax=axs[2,0]);\r\n",
    "l2, r2 = axs[2, 0].get_xlim()\r\n",
    "axs[2, 1].set_xlim(left=l2, right=r2)\r\n",
    "sns.boxplot(x=\"value\", y=\"variable\", data=pd.melt(boxp3), showmeans=True, palette=\"Set2\", ax=axs[2,1]);\r\n",
    "\r\n",
    "# just for fun let's add a plot of the actual data to the last row\r\n",
    "sns.stripplot(x=\"value\", y=\"variable\", data=pd.melt(boxp3), marker=\"o\", alpha=0.3, color=\"black\", ax=axs[2,0]);\r\n",
    "sns.stripplot(x=\"value\", y=\"variable\", data=pd.melt(boxp3), marker=\"o\", alpha=0.3, color=\"black\", ax=axs[2,1]);"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Everything I have read says the value of boxplots is in making quick comparisons. Not show much in looking at features in a single sample. Histograms and such are likely better for the latter case. Though the violin plot sort of combines some of the features of both.\r\n",
    "\r\n",
    "So, let's have look at some possible comparisons. For example BMI and/or BP by sex; as we did with histograms above."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "fig, axs = plt.subplots(2, 2, figsize=(10, 10))\r\n",
    "cnt_x1 = db_data[db_data['SEX'] == 1].SEX.count()\r\n",
    "cnt_x2 = db_data[db_data['SEX'] == 2].SEX.count()\r\n",
    "print(f\"sex 1: {cnt_x1}; sex 2: {cnt_x2}\")\r\n",
    "#sns.boxplot(x='BMI', y='Y', data=boxp3, notch=True, hue=\"SEX\")\r\n",
    "#sns.boxplot(x=\"value\", y=\"variable\", data=pd.melt(boxp2), showmeans=True, orient=\"h\", palette=\"Set3\", ax=axs[1,1]);\r\n",
    "sns.boxplot(x=\"SEX\", y=\"BMI\", data=db_data, palette=\"Set3\", notch=True, ax=axs[0, 0]);\r\n",
    "sns.boxplot(x=\"SEX\", y=\"BP\", data=db_data, palette=\"Set3\", notch=True, ax=axs[0, 1]);\r\n",
    "sns.boxplot(x=\"SEX\", y=\"AGE\", data=db_data, palette=\"Set3\", notch=True, ax=axs[1,0]);\r\n",
    "sns.boxplot(x=\"SEX\", y=\"S4\", data=db_data, palette=\"Set3\", notch=True, ax=axs[1,1]);\r\n",
    "axs[1, 1].set_ylabel(\"TSH (tch)\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Ok! That is possibly more interesting.\r\n",
    "\r\n",
    "The median BMI for both sexes is similar. And, as the notches overlap, we definitely can not say that the true medians are different. But sex 1 appears to have more variablility in BMI than sex 2. With the sex 2 distribution perhaps a bit more skewed.\r\n",
    "\r\n",
    "With respect to blood pressure, sex 2 has a higher median pressure. And, as the notches do not overlap, we can conclude, with 95% confidence, that the true medians do differ. The variability between the sexes seems similar. Perhaps a bit more skew for sex 1.\r\n",
    "\r\n",
    "Sex 2 has a higher median age. And the true median age for the sexes would appear to be different. Variability once again similar.\r\n",
    "\r\n",
    "TSH certainly looks stranger. Would appear to very definitely be some skewness in the distribution for both sexes. Again variability similar. But the mean TSH for sex 2 definitely higher than for sex 1. Looks like almost 75% of sex 2 have a TSH higher tham 75% of sex 1.\r\n",
    "\r\n",
    "Not sure this would help with any attempt to model the data, but&#8212;"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# Let's look at some bivariate visualization\r\n",
    "#plt.figure();\r\n",
    "g = sns.relplot(x=\"BMI\", y=\"Y\", hue=\"SEX\", data=db_data, palette=\"Set2\");\r\n",
    "#plt.ylabel(\"Progression\")\r\n",
    "g.set_ylabels(\"Progression\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "g = sns.lmplot(x=\"BMI\", y=\"Y\", hue=\"SEX\", data=db_data, palette=\"Set2\");\r\n",
    "#plt.ylabel(\"Progression\")\r\n",
    "g.set_ylabels(\"Progression\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "g = sns.lmplot(x=\"BMI\", y=\"Y\", col=\"SEX\", data=db_data, palette=\"Set2\");\r\n",
    "#plt.ylabel(\"Progression\")\r\n",
    "g.set_ylabels(\"Progression\");"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# also look at seaborn tips"
   ],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "orig_nbformat": 4,
  "language_info": {
   "name": "python",
   "version": "3.9.2",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.2 64-bit ('ds-3.9': conda)"
  },
  "interpreter": {
   "hash": "a27d3f2bf68df5402465348834a2195030d3fc5bfc8e594e2a17c8c7e2447c85"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}