{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "078900c6-d71d-44f0-9e76-cc5ce66f0ca9",
   "metadata": {},
   "source": [
    "# Logistic regression with uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0de7f02-942b-42ff-874c-3c0e2ac1e0dc",
   "metadata": {},
   "source": [
    "## Data & common "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "be75450c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from fitter import Fitter, get_common_distributions, get_distributions\n",
    "from scipy import io\n",
    "\n",
    "import scipy.stats as ss\n",
    "\n",
    "from sklearn import linear_model\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3373f87f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#data_path = \"/home/marija/maki/\"\n",
    "data_path = \"../data/\"\n",
    "suv_filename = data_path + \"suv_percentilesSLOthenUWM.mat\"\n",
    "flags_filename=  data_path + \"flags_combined.mat\"\n",
    "normal_range_filename =  data_path + \"normal_range.mat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "811d990d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Loading all data\n",
    "\"\"\"\n",
    "suv_dict = io.loadmat(suv_filename)\n",
    "flags_dict = io.loadmat(flags_filename)\n",
    "normal_range_dict = io.loadmat(normal_range_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b47b0289-1378-40a8-88e2-acfed4c4cd69",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4",
   "metadata": {},
   "source": [
    "### General ideas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ae5d1e06-5122-400d-a9ab-7f04867f8cdb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 17, 100)\n",
      "(58,)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Extracting only data associated to organ of interest, you need to choose data and column of flags that shoud be taken depended on data \n",
    "\"\"\"\n",
    "# in suv values can be (choose one of them): 'lung_SUVperc_COMBINED' , 'bowel_SUVperc_COMBINED' or 'thyroid_SUVperc_COMBINED'\n",
    "# in flags choose the column depended on the data of interest (1 for bowel, 3 for lung, 5 for thyroid)\n",
    "\n",
    "suv = suv_dict['lung_SUVperc_COMBINED'][0:58,:,:]\n",
    "flags = flags_dict['flags'][0:58,3]\n",
    "\n",
    "print(suv.shape)\n",
    "print(flags.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a70c64a3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      "[1 1 1 1 1]\n"
     ]
    }
   ],
   "source": [
    "print(flags[flags==0])\n",
    "print(flags[flags==1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9b894e29",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(58, 1)\n"
     ]
    }
   ],
   "source": [
    "p = 94  # choose percentile of interest (for lung and bowel 95, for thyroid 75), numerating here is from 0 so take the one less\n",
    "X = np.nanmax(suv[:,:,p], axis=1).reshape(-1,1)\n",
    "\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6deba69b-9c9c-47b6-b02d-f40998729aa1",
   "metadata": {},
   "source": [
    "Logistic regression documentation:\n",
    "* https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression\n",
    "* https://stats.stackexchange.com/questions/186830/what-is-scikit-learns-logisticregression-minimizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "239ae789",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w0: [-11.68397793]\n",
      "w1: [[5.68702693]]\n"
     ]
    }
   ],
   "source": [
    "logr = linear_model.LogisticRegression(penalty='none')\n",
    "clf = logr.fit(X, flags)\n",
    "\n",
    "# underlaying linear model coefficients\n",
    "print(\"w0:\", logr.intercept_)\n",
    "print(\"w1:\", logr.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "54f463d0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb6b21c2560>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 19096 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 6, 100)\n",
    "y = logr.predict_proba(x.reshape(-1,1))\n",
    "\n",
    "plt.title(\"Prediction probability of AE depending on maxSUV value\")\n",
    "plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "plt.ylabel(\"P(AE|X=x)\")\n",
    "\n",
    "plt.scatter(X[flags==0], flags[flags==0],label=\"NC\")  \n",
    "plt.scatter(X[flags==1], flags[flags==1],label=\"AE\")\n",
    "\n",
    "plt.plot(x, 1/(1 + np.exp(-logr.intercept_[0] - logr.coef_[0,0]*x)), label=\"logistic reg\")\n",
    "plt.legend(loc='center right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2c16aad9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covariance matrix: \n",
      " [[12.32329923 -6.71305958]\n",
      " [-6.71305958  3.81886093]]\n",
      "Standard errors:  [3.5104557  1.95419061]\n",
      "logitParams: [-11.68397793   5.68702693]\n"
     ]
    }
   ],
   "source": [
    "#Ref: \n",
    "# https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients\n",
    "\n",
    "# Calculate matrix of predicted class probabilities.\n",
    "predProbs = logr.predict_proba(X) #has two col. p, (1-p)\n",
    "\n",
    "# Design matrix -- add column of 1's at the beginning of your X_train matrix\n",
    "X_design = np.hstack([np.ones((X.shape[0], 1)), X]) \n",
    "\n",
    "# Initiate matrix of 0's, fill diagonal with each predicted observation's variance\n",
    "V = np.diagflat(np.product(predProbs, axis=1)) #dig.matrix where each element is p*(1-p)\n",
    "\n",
    "# Covariance matrix C_params = (X^T V X)^-1, params = (b0,b1) \n",
    "covLogit = np.linalg.inv(np.dot(np.dot(X_design.T, V), X_design))\n",
    "print(\"Covariance matrix: \\n\", covLogit)\n",
    "\n",
    "# Standard errors\n",
    "print(\"Standard errors: \", np.sqrt(np.diag(covLogit)))\n",
    "\n",
    "# Wald statistic (coefficient / s.e.) ^ 2\n",
    "logitParams = np.insert(logr.coef_, 0, logr.intercept_)\n",
    "print(\"logitParams:\", logitParams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2f0ac198",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Logistic regresion modeling with estimated uncertainties')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 24336 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Direct approach to variance in logit regression:\n",
    "#\n",
    "#   f(x|p*)    p ~ N(p*,C_p)\n",
    "#\n",
    "\n",
    "f = lambda x,p: 1/(1 + np.exp(-p[0] - p[1]*x))\n",
    "fopt = lambda x: 1/(1 + np.exp(-logitParams[0] - logitParams[1]*x))\n",
    "\n",
    "#generate parameters and points in x axis\n",
    "ps = np.random.multivariate_normal(logitParams, covLogit, 1000)\n",
    "xs = np.linspace(0, 6, 100)\n",
    "\n",
    "# generate values in y axis\n",
    "yopts = fopt(xs)\n",
    "ys = np.array([np.fromiter(map(lambda p: f(x,p), ps),float) for x in xs])\n",
    "\n",
    "# doing statistics on values\n",
    "mys = np.median(ys,axis=1)\n",
    "alpha=0.05\n",
    "yerr_min = np.quantile(ys, alpha/2, axis=1)\n",
    "yerr_max = np.quantile(ys, 1-alpha/2, axis=1)\n",
    "\n",
    "plt.xlabel(\"x = maxSUV\")\n",
    "plt.ylabel(\"P(AE|X=x)\")\n",
    "plt.scatter(X[flags==0], flags[flags==0],label=\"NC\")  \n",
    "plt.scatter(X[flags==1], flags[flags==1],label=\"AE\")\n",
    "plt.plot(xs, yopts, label=\"fit logit model\",color=\"red\")\n",
    "#plt.errorbar(xs, mys, yerr = [mys-yerr_min, yerr_max-mys], label=\"direct variance\",fmt=' ')\n",
    "plt.fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"confidence interval\")\n",
    "\n",
    "plt.legend(loc='center right')\n",
    "plt.title(\"Logistic regresion modeling with estimated uncertainties\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a1a899c-ac8a-44e8-bad6-04262e429913",
   "metadata": {},
   "source": [
    "### Distribution of parameters via BS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "712f33c6-b772-4c90-b2e0-ad9336fc8d83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading data for lungs\n",
    "suv = suv_dict['lung_SUVperc_COMBINED'][0:58,:,:]\n",
    "flags = flags_dict['flags'][0:58,3]\n",
    "\n",
    "per = 94  # choose percentile of interest (for lung and bowel 95, for thyroid 75), numerating here is from 0 so take the one less\n",
    "X = np.nanmax(suv[:,:,per], axis=1).reshape(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8d8c5f8a-f065-4980-a313-786e0df9b2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss = lambda x, mu, var: np.exp(-(x - mu)**2/(2*var))/np.sqrt(2*np.pi*var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d2ff06ca-e6f3-4442-9bce-29ee929c0071",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()\n",
    "\n",
    "def get_params(idx = None):\n",
    "    if idx is None:\n",
    "        idx = rng.choice(n, n)\n",
    "        while (flags[idx] == 0).all() or (flags[idx] == 1).all():\n",
    "            idx = rng.choice(n, n)\n",
    "            \n",
    "    clf = linear_model.LogisticRegression(penalty='none').fit(X[idx], flags[idx])\n",
    "    return np.insert(clf.coef_, 0, clf.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "69792c6a-18d5-43c1-b1f9-a0744beb1ebb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# original p0\n",
    "n = len(X)\n",
    "p0 = get_params(np.arange(n))\n",
    "\n",
    "# bootstrapping\n",
    "ps = np.array([get_params() for _ in range(10**7)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a97debd2-00cf-4ebc-9b50-da847daf6d63",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27576 characters]\n",
      "text/plain": [
       "<Figure size 600x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting distribution of parameters\n",
    "sns.set(rc={\"figure.figsize\":(6, 4), \"font.size\": 14})\n",
    "\n",
    "ax = sns.jointplot(x = ps[:,0], y = ps[:,1], kind ='kde', fill=True)\n",
    "ax.set_axis_labels(xlabel= f\"$\\\\beta_0$\", ylabel=f\"$\\\\beta_1$\")\n",
    "\n",
    "plt.text(0.04, 0.9, \"(a)\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=ax.figure.transFigure)\n",
    "plt.savefig(\"./figs/logit_boots.pdf\", bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4f08a287-fdcb-47cc-a2fd-55ddba066ce8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size_selectd: 6298713\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 48808 characters]\n",
      "text/plain": [
       "<Figure size 600x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# select main peak and check distribution\n",
    "ps_sel = ps[(ps[:,0] > -25) & ( ps[:,0] < 0) & (ps[:,1]> 0) & ( ps[:,1] < 10)];\n",
    "print(\"size_selectd:\", len(ps_sel))\n",
    "\n",
    "# plotting distribution of parameters\n",
    "sns.set(rc={\"figure.figsize\":(6, 4), \"font.size\":14})\n",
    "\n",
    "ax = sns.jointplot(x = ps_sel[:,0], y = ps_sel[:,1], kind ='kde', fill=True)\n",
    "ax.set_axis_labels(xlabel= f\"$\\\\beta_0$\", ylabel=f\"$\\\\beta_1$\")\n",
    "\n",
    "plt.text(0.04, 0.9, \"(b)\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=ax.figure.transFigure)\n",
    "plt.plot([p0[0]], [p0[1]], marker=\"o\", markersize=10, markeredgecolor=\"red\", markerfacecolor=\"red\")\n",
    "\n",
    "x = np.linspace(-25,0,100)\n",
    "ax.ax_marg_x.plot(x, gauss(x, logitParams[0], covLogit[0][0]), \"r\" )\n",
    "\n",
    "y = np.linspace(0, 10,100)\n",
    "ax.ax_marg_y.plot(gauss(y, logitParams[1], covLogit[1][1]), y, \"r\" )\n",
    "\n",
    "plt.savefig(\"./figs/logit_boots_sel.pdf\", bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "668f5b02-e351-4dde-9697-af36b49588f2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLE: beta:\n",
      " [-11.68397793   5.68702693]\n",
      "MLE: cov(beta):\n",
      " [[12.32329923 -6.71305958]\n",
      " [-6.71305958  3.81886093]]\n"
     ]
    }
   ],
   "source": [
    "# MLE: optimal beta and covariance\n",
    "print(\"MLE: beta:\\n\", logitParams)\n",
    "print(\"MLE: cov(beta):\\n\", covLogit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "308b94a4-7a95-47c8-b0c5-4a80bc05b3df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample_mean(beta):\n",
      " [-10.82981446   5.00156176]\n",
      "sample_cov(beta):\n",
      " [[ 4.5526891  -2.65783521]\n",
      " [-2.65783521  1.70060809]]\n"
     ]
    }
   ],
   "source": [
    "# sample mean and covariance of selected main peak\n",
    "smean = np.mean(ps_sel, axis=0)\n",
    "print(\"sample_mean(beta):\\n\", smean)\n",
    "db = ps_sel - smean\n",
    "scov = np.dot(db.T, db)/(len(db)-1)\n",
    "print(\"sample_cov(beta):\\n\", scov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "cc3d41a7-518f-4804-832a-090cd1501d35",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35356 characters]\n",
      "text/plain": [
       "<Figure size 600x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# zooming specific area and marking MLE solution:\n",
    "sns.set(rc={\"figure.figsize\":(6, 4), \"font.size\":14})\n",
    "\n",
    "ax = sns.jointplot(x = ps[:,0], y = ps[:,1], kind ='kde', fill=True)\n",
    "ax.set_axis_labels(xlabel= f\"$\\\\beta_0$\", ylabel=f\"$\\\\beta_1$\")\n",
    "\n",
    "x = np.linspace(-25,10,100)\n",
    "ax.ax_marg_x.plot(x, gauss(x, logitParams[0], covLogit[0][0]), \"r\" )\n",
    "\n",
    "y = np.linspace(-20,40,100)\n",
    "ax.ax_marg_y.plot(gauss(y, logitParams[1], covLogit[1][1]), y, \"r\" )\n",
    "\n",
    "plt.text(0.04, 0.9, \"(b)\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=ax.figure.transFigure)\n",
    "plt.plot([p0[0]], [p0[1]], marker=\"o\", markersize=10, markeredgecolor=\"red\", markerfacecolor=\"red\")\n",
    "\n",
    "plt.xlim(-80,30)\n",
    "plt.ylim(-20,40)\n",
    "plt.savefig(\"./figs/logit_boots_zoom.pdf\", bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4c2f7492-6caa-4f05-8a2f-a5a1eac28d30",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact: opt(p)= [-11.68397793   5.68702693]\n",
      "exact: covlogit= [[12.32329923 -6.71305958]\n",
      " [-6.71305958  3.81886093]]\n",
      "mean(p)= [-45.79078382  25.75168304]\n",
      "cov(p)= [[ 5276.30572959 -3232.9241108 ]\n",
      " [-3232.9241108   1991.14599825]]\n"
     ]
    }
   ],
   "source": [
    "# calculate empirical average and variance-covariance matrix\n",
    "pmean = np.mean(ps, axis=0)\n",
    "dp = ps - pmean\n",
    "pcov = np.dot(dp.T, dp)/(len(dp)-1)\n",
    "\n",
    "print(\"exact: opt(p)=\", p0)\n",
    "print(\"exact: covlogit=\",covLogit)\n",
    "\n",
    "print(\"mean(p)=\", pmean)\n",
    "print(\"cov(p)=\", pcov)\n",
    "\n",
    "# Conclusion: due to multiple peaks the empirical values are very far from MLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9bd0c944-1793-42f5-8a78-8a38e339e8f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Install:\n",
    "#!pip install pingouin\n",
    "#\n",
    "# Ref: \n",
    "#  * https://www.statology.org/multivariate-normality-test-python/\n",
    "#  * https://pingouin-stats.org/build/html/generated/pingouin.multivariate_normality.html\n",
    "# https://www.sfu.ca/sasdoc/sashtml/ets/chap14/sect38.htm\n",
    "\n",
    "from pingouin import multivariate_normality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5160ac22-29e7-4c06-8b60-1214231fbc43",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HZResults(hz=0.6501826885176032, pval=0.6851557047523155, normal=True)\n",
      "HZResults(hz=10.665573669445136, pval=7.185502334798263e-24, normal=False)\n"
     ]
    }
   ],
   "source": [
    "# testing multivariate_normality\n",
    "want_to_test = False\n",
    "\n",
    "if want_to_test:\n",
    "    # test cases\n",
    "    test1 = multivariate_normality(rng.multivariate_normal([1,2],[[2,0.5],[0.5,1]], size = 1000), alpha=.05)\n",
    "    \n",
    "    data = rng.multivariate_normal([3,5],[[1,0.5],[0.5,1]], size = 1000)\n",
    "    test2 = multivariate_normality(np.log(np.abs(data)), alpha=.05)\n",
    "    \n",
    "    print(test1)\n",
    "    print(test2)\n",
    "\n",
    "if want_to_test:\n",
    "    # testing our paramaters\n",
    "    our_case = multivariate_normality(ps, alpha=.05)\n",
    "    print(our_case)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98ad2d95-d89b-4f18-a2f7-c7fc17473135",
   "metadata": {},
   "source": [
    "Conclusion: Statistically looking, the distribution of parameters is far from multivariate Gaussian for these sizes of samples, but is captures essential properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e065c5d1-7dd5-4098-a095-5a283a8d408c",
   "metadata": {},
   "source": [
    "## Final results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8305144d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "figsdir = \"./figs/\"\n",
    "\n",
    "def gen_plots(pert, organ, prefix, alpha, title, bootstrapping = False):\n",
    "    \n",
    "    suv = suv_dict[organ + '_SUVperc_COMBINED'][0:58,:,:]\n",
    "  \n",
    "    if organ == \"lung\":\n",
    "        idx = 3\n",
    "    elif organ == \"bowel\":\n",
    "        idx = 1\n",
    "    elif organ == \"thyroid\":\n",
    "        idx = 5\n",
    "        \n",
    "    flags = flags_dict['flags'][0:58, idx]\n",
    "    \n",
    "    # choose percentile of interest (for lung and bowel 95, for thyroid 75), \n",
    "    # numerating here is from 0 so take the one less\n",
    "    X = np.nanmax(suv[:, :, pert - 1], axis=1).reshape(-1,1)\n",
    "    \n",
    "    # fitting\n",
    "    clf = linear_model.LogisticRegression(penalty='none').fit(X, flags)\n",
    "\n",
    "    # Store fitted parameters in array\n",
    "    logitParams = np.insert(clf.coef_, 0, clf.intercept_)\n",
    "    \n",
    "    # logit function\n",
    "    f = lambda x, p: 1/(1 + np.exp(-p[0] - p[1]*x))\n",
    "\n",
    "    Nsamples = 10**5\n",
    "    \n",
    "    if bootstrapping:\n",
    "        # obtaining parameters via bootstrapping method\n",
    "        rng = np.random.default_rng()\n",
    "\n",
    "        def BS_sample():\n",
    "            idx = rng.integers(len(X), size=len(X))\n",
    "            st = (flags[idx] == 0)\n",
    "            \n",
    "            while st.all() or (~st).all(): \n",
    "                idx = rng.integers(len(X), size=len(X))\n",
    "                st = (flags[idx] == 0)\n",
    "                \n",
    "            tmp = linear_model.LogisticRegression(penalty='none').fit(X[idx], flags[idx])\n",
    "            \n",
    "            return np.insert(tmp.coef_, 0, tmp.intercept_)\n",
    "        \n",
    "        ps = np.array([BS_sample() for _ in range(Nsamples)])\n",
    "\n",
    "    else:\n",
    "        # Calculate matrix of predicted class probabilities.\n",
    "        predProbs = clf.predict_proba(X) #has two col. p, (1-p)\n",
    "        \n",
    "        # Design matrix -- add column of 1's at the beginning of your X_train matrix\n",
    "        X_design = np.hstack([np.ones((X.shape[0], 1)), X]) \n",
    "    \n",
    "        # Initiate matrix of 0's, fill diagonal with each predicted observation's variance\n",
    "        V = np.diagflat(np.prod(predProbs, axis=1)) #dig.matrix where each element is p*(1-p)\n",
    "    \n",
    "        # Covariance matrix C_params = (X^T V X)^-1, params = (b0,b1) \n",
    "        covLogit = np.linalg.inv(np.dot(np.dot(X_design.T, V), X_design))\n",
    "        \n",
    "        #generate parameters and points in x axis\n",
    "        ps = np.random.multivariate_normal(logitParams, covLogit, Nsamples)\n",
    "        \n",
    "    # get for plotting purpose meaningful xmax\n",
    "    # f = 1- eps, np.exp(-p[0] - p[1]*x) ~ eps => -log(eps) = p[0] + p[1]*x\n",
    "    xmax = max(X.max(), (-np.log(1e-3)-logitParams[0])/logitParams[1])\n",
    "    \n",
    "    xs = np.linspace(0, xmax, 200)\n",
    "    \n",
    "    # generate values in y axis\n",
    "    yopts = f(xs, logitParams)\n",
    "    ys = np.array([np.fromiter(map(lambda p: f(x,p), ps),float) for x in xs])\n",
    "\n",
    "    # calculate median and quantiles\n",
    "    mys = np.median(ys, axis=1)\n",
    "    yerr_min = np.quantile(ys, alpha/2, axis=1)\n",
    "    yerr_max = np.quantile(ys, 1-alpha/2, axis=1)\n",
    "\n",
    "    # generate plots\n",
    "    plt.clf()\n",
    "    fs=14\n",
    "    params = {'font.size': fs, 'axes.labelsize': fs, 'axes.titlesize':fs, 'legend.fontsize': fs, 'xtick.labelsize': fs, 'ytick.labelsize': fs}\n",
    "    plt.rcParams.update(params)\n",
    "\n",
    "    fig = plt.figure(figsize = (7, 5))\n",
    "    \n",
    "    #plt.xlabel(r\"$x = max_{visit} SUV(visit, p)$\")\n",
    "    #plt.ylabel(\"P(YAE|X=x)\")\n",
    "    \n",
    "    plt.xlabel(r\"$x$\")\n",
    "    plt.ylabel(r\"$P(Y = AE|x)$\")\n",
    "   \n",
    "    plt.scatter(X[flags==0], flags[flags==0], label=\"NC\")  \n",
    "    plt.scatter(X[flags==1], flags[flags==1], label=\"AE\")\n",
    "    \n",
    "    plt.plot(xs, yopts, label=\"optimal\", color=\"red\")\n",
    "    plt.plot(xs, mys, label=\"median\", color=\"blue\")\n",
    "    \n",
    "    plt.fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"confidence band\")\n",
    "\n",
    "    #plt.title(title+\":\" + organ)\n",
    "    #plt.text()\n",
    "    plt.text(0.1, 0.93, \"(\"+organ+\")\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=fig.transFigure)\n",
    "    \n",
    "    plt.legend(loc = 'lower right', fancybox=True, fontsize=fs, framealpha=0.5)\n",
    "    plt.tight_layout()\n",
    "\n",
    "    if not os.path.exists(figsdir): os.makedirs(figsdir)\n",
    "        \n",
    "    plt.savefig(figsdir + prefix + \"_\" + organ + \".pdf\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def gen_allplots(perts, organs, prefix, alpha, bootstrapping = False):\n",
    "\n",
    "    fs = 14\n",
    "    params = {'font.size': fs, 'axes.labelsize': fs, 'axes.titlesize':fs, 'legend.fontsize': fs, 'xtick.labelsize': fs, 'ytick.labelsize': fs}\n",
    "    plt.rcParams.update(params)\n",
    "    \n",
    "    fig, axs = plt.subplots(nrows=1, ncols=len(organs), sharey = True, figsize=(6*len(organs), 4))\n",
    "    \n",
    "    res_fit_pars = []\n",
    "    \n",
    "    for i, organ in enumerate(organs):\n",
    "        suv = suv_dict[organ + '_SUVperc_COMBINED'][0:58,:,:]\n",
    "      \n",
    "        if organ == \"lung\":\n",
    "            idx = 3\n",
    "        elif organ == \"bowel\":\n",
    "            idx = 1\n",
    "        elif organ == \"thyroid\":\n",
    "            idx = 5\n",
    "            \n",
    "        flags = flags_dict['flags'][0:58, idx]\n",
    "        pert = perts[i]\n",
    "        \n",
    "        # choose percentile of interest (for lung and bowel 95, for thyroid 75), \n",
    "        # numerating here is from 0 so take the one less\n",
    "        X = np.nanmax(suv[:, :, pert - 1], axis=1).reshape(-1,1)\n",
    "        \n",
    "        # fitting\n",
    "        clf = linear_model.LogisticRegression(penalty='none').fit(X, flags)\n",
    "    \n",
    "        # Store fitted parameters in array\n",
    "        logitParams = np.insert(clf.coef_, 0, clf.intercept_)\n",
    "\n",
    "        # store fit result in res\n",
    "        res_fit_pars.append((organ, logitParams[0],  logitParams[1]))\n",
    "        \n",
    "        # logit function\n",
    "        f = lambda x, p: 1/(1 + np.exp(-p[0] - p[1]*x))\n",
    "    \n",
    "        Nsamples = 10**5\n",
    "        \n",
    "        if bootstrapping:\n",
    "            # obtaining parameters via bootstrapping method\n",
    "            rng = np.random.default_rng()\n",
    "    \n",
    "            def BS_sample():\n",
    "                idx = rng.integers(len(X), size=len(X))\n",
    "                st = (flags[idx] == 0)\n",
    "                \n",
    "                while st.all() or (~st).all(): \n",
    "                    idx = rng.integers(len(X), size=len(X))\n",
    "                    st = (flags[idx] == 0)\n",
    "                    \n",
    "                tmp = linear_model.LogisticRegression(penalty='none').fit(X[idx], flags[idx])\n",
    "                \n",
    "                return np.insert(tmp.coef_, 0, tmp.intercept_)\n",
    "            \n",
    "            ps = np.array([BS_sample() for _ in range(Nsamples)])\n",
    "    \n",
    "        else:\n",
    "            # Calculate matrix of predicted class probabilities.\n",
    "            predProbs = clf.predict_proba(X) #has two col. p, (1-p)\n",
    "            \n",
    "            # Design matrix -- add column of 1's at the beginning of your X_train matrix\n",
    "            X_design = np.hstack([np.ones((X.shape[0], 1)), X]) \n",
    "        \n",
    "            # Initiate matrix of 0's, fill diagonal with each predicted observation's variance\n",
    "            V = np.diagflat(np.prod(predProbs, axis=1)) #dig.matrix where each element is p*(1-p)\n",
    "        \n",
    "            # Covariance matrix C_params = (X^T V X)^-1, params = (b0,b1) \n",
    "            covLogit = np.linalg.inv(np.dot(np.dot(X_design.T, V), X_design))\n",
    "            \n",
    "            #generate parameters and points in x axis\n",
    "            ps = np.random.multivariate_normal(logitParams, covLogit, Nsamples)\n",
    "            \n",
    "        # get for plotting purpose meaningful xmax\n",
    "        # f = 1- eps, np.exp(-p[0] - p[1]*x) ~ eps => -log(eps) = p[0] + p[1]*x\n",
    "        xmax = max(X.max(), (-np.log(1e-3)-logitParams[0])/logitParams[1])\n",
    "        \n",
    "        xs = np.linspace(0, xmax, 200)\n",
    "        \n",
    "        # generate values in y axis\n",
    "        yopts = f(xs, logitParams)\n",
    "        ys = np.array([np.fromiter(map(lambda p: f(x,p), ps),float) for x in xs])\n",
    "    \n",
    "        # calculate median and quantiles\n",
    "        mys = np.median(ys, axis=1)\n",
    "        yerr_min = np.quantile(ys, alpha/2, axis=1)\n",
    "        yerr_max = np.quantile(ys, 1-alpha/2, axis=1)\n",
    "\n",
    "        axs[i].set_xlabel(r\"$x$\")\n",
    "        if i == 0: axs[i].set_ylabel(r\"$P(Y = AE|x)$\")\n",
    "    \n",
    "        axs[i].scatter(X[flags==0], flags[flags==0], label=\"NC\")  \n",
    "        axs[i].scatter(X[flags==1], flags[flags==1], label=\"AE\")\n",
    "        \n",
    "        axs[i].plot(xs, yopts, label=\"optimal\", color=\"red\")\n",
    "        axs[i].plot(xs, mys, label=\"median\", color=\"blue\")\n",
    "        \n",
    "        axs[i].fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"confidence band\")\n",
    "        axs[i].text(0.1, 0.90, \"(\"+organ+\")\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=axs[i].transAxes)\n",
    "\n",
    "    #plt.title(title+\":\" + organ)\n",
    "    #plt.text()\n",
    "    #plt.text(0.1, 0.93, \"(\"+organ+\")\", fontsize=16, horizontalalignment='center', verticalalignment='center', transform=fig.transFigure)\n",
    "\n",
    "    if not os.path.exists(figsdir): os.makedirs(figsdir)\n",
    "    \n",
    "    loc = 'lower center'\n",
    "    ncol = 5\n",
    "    handles, labels = axs[-1].get_legend_handles_labels()\n",
    "    fig.legend(handles=handles, labels=labels, loc=loc, ncol=ncol, bbox_to_anchor=(0.5, -0.1))    \n",
    "    \n",
    "    plt.subplots_adjust(wspace=0, hspace=0)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(figsdir + prefix + \".pdf\", bbox_inches='tight')\n",
    "    plt.show()\n",
    "\n",
    "    # print results\n",
    "    print(pd.DataFrame(res_fit_pars, columns=[\"organ\", \"beta0\", \"beta1\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb7dda9a-2deb-44ac-823d-66773165be59",
   "metadata": {},
   "source": [
    "### Multivaiate normal approx of parameter pdf -- \"parametric BS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b3a0f89-bad4-4a39-81b9-1b9ac8ca3f75",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "title=\"Logistic regresion with estimated uncertainties\"\n",
    "gen_plots(95, \"lung\", \"logit_mvn\", 0.05, title)\n",
    "gen_plots(95, \"bowel\", \"logit_mvn\", 0.05, title)\n",
    "gen_plots(75, \"thyroid\", \"logit_mvn\", 0.05, title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d5a1c0-6b4f-40c1-b866-4cd0baa05159",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.rcParams.keys() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f116b595-6fc5-47cb-80be-21b50b3422e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 60580 characters]\n",
      "text/plain": [
       "<Figure size 1296x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     organ      beta0     beta1\n",
      "0     lung -11.683978  5.687027\n",
      "1    bowel  -3.474206  0.337083\n",
      "2  thyroid  -4.387076  1.163232\n"
     ]
    }
   ],
   "source": [
    "gen_allplots([95, 95, 75], [\"lung\", \"bowel\", \"thyroid\"], \"logit_mvn\", 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1e9c470-9a7d-4490-9912-be9d24a31bab",
   "metadata": {},
   "source": [
    "### Non-parametric Bootstrapping (BS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c461077-18cd-4200-88fb-f10976449eb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "title=\"Logistic regresion with estimated uncertainties\"\n",
    "gen_plots(95, \"lung\", \"logit_bs\", 0.05, title, bootstrapping = True)\n",
    "gen_plots(95, \"bowel\", \"logit_bs\", 0.05, title, bootstrapping = True)\n",
    "gen_plots(75, \"thyroid\", \"logit_bs\", 0.05, title, bootstrapping = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20503b2a-51c4-4ad1-8cff-d222d734e198",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
