{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "078900c6-d71d-44f0-9e76-cc5ce66f0ca9",
   "metadata": {
    "id": "078900c6-d71d-44f0-9e76-cc5ce66f0ca9"
   },
   "source": [
    "# Logistic regression with uncertainties\n",
    "\n",
    "Description:\n",
    "  * Two parameter model\n",
    "  * To obtain confidence intervals of model values, two types of bootstrapping are performed:\n",
    "    * a bit constrained non-parametric bootstrapping\n",
    "    * parametric bootstrapping of parameters around highest peak obtained via non-parametric bootstrapping"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0de7f02-942b-42ff-874c-3c0e2ac1e0dc",
   "metadata": {
    "id": "e0de7f02-942b-42ff-874c-3c0e2ac1e0dc"
   },
   "source": [
    "## Data & common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "40c9c487-ff42-4463-b105-4d7bf75d8f6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "choice_setup = 1  # 0 for MD, 1 for MH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "56ecf5e9-308d-411b-be22-c5a1ac001717",
   "metadata": {},
   "outputs": [],
   "source": [
    "if choice_setup == 0:\n",
    "    from google.colab import drive\n",
    "    drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "be75450c",
   "metadata": {
    "id": "be75450c",
    "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",
    "import scipy\n",
    "from scipy import io\n",
    "import scipy.stats as ss\n",
    "\n",
    "from sklearn import linear_model\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3373f87f",
   "metadata": {
    "id": "3373f87f",
    "tags": []
   },
   "outputs": [],
   "source": [
    "topic = \"logit\"\n",
    "\n",
    "if choice_setup == 0:\n",
    "    #data_path = \"/home/marija/maki/\"\n",
    "    data_path = \"drive/MyDrive/maki_slovenija/data/\"\n",
    "    figsdir = f\"drive/MyDrive/maki_slovenija/figs_{topic}/\"\n",
    "    figsext = \".eps\"  # figure extension\n",
    "elif choice_setup == 1:\n",
    "    data_path = \"../data\"\n",
    "    figsdir = f\"./figs/{topic}\"\n",
    "    figsext = \".pdf\"  # figure extension\n",
    "    \n",
    "if not os.path.exists(figsdir): os.makedirs(figsdir)\n",
    "        \n",
    "suv_filename = os.path.join(data_path, \"suv_percentilesSLOthenUWM.mat\")\n",
    "flags_filename = os.path.join(data_path,  \"flags_combined.mat\")\n",
    "normal_range_filename = os.path.join(data_path, \"normal_range.mat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "811d990d",
   "metadata": {
    "id": "811d990d",
    "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": {
    "id": "b47b0289-1378-40a8-88e2-acfed4c4cd69"
   },
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4",
   "metadata": {
    "id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4"
   },
   "source": [
    "### General ideas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ae5d1e06-5122-400d-a9ab-7f04867f8cdb",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 9,
     "status": "ok",
     "timestamp": 1733838433084,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "ae5d1e06-5122-400d-a9ab-7f04867f8cdb",
    "outputId": "a55e5243-c612-4fa3-9da9-662a064b45e5"
   },
   "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": 7,
   "id": "a70c64a3",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 8,
     "status": "ok",
     "timestamp": 1733838433084,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "a70c64a3",
    "outputId": "a4cc1e1b-e752-45b5-a76d-43d3f664b159",
    "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": 8,
   "id": "9b894e29",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1733838433084,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "9b894e29",
    "outputId": "1797f635-0fc6-46af-8a3f-6935f385383b",
    "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": {
    "id": "6deba69b-9c9c-47b6-b02d-f40998729aa1"
   },
   "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": "markdown",
   "id": "FhyIAER1E7PZ",
   "metadata": {
    "id": "FhyIAER1E7PZ"
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "239ae789",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 492,
     "status": "ok",
     "timestamp": 1733842169917,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "239ae789",
    "outputId": "ed8c1661-233c-47a7-b2a1-136677e6fe36",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w0: [-5.33495652]\n",
      "w1: [[1.90942694]]\n"
     ]
    }
   ],
   "source": [
    "#ovo nam ne odgovara jer nama treba regularizacija da bi nam parametri bili ograniceni\n",
    "#logr = linear_model.LogisticRegression(penalty='none')\n",
    "#clf = logr.fit(X, flags)\n",
    "clf = linear_model.LogisticRegression(random_state=0).fit(X, flags)\n",
    "\n",
    "# underlaying linear model coefficients\n",
    "#print(\"w0:\", logr.intercept_)\n",
    "#print(\"w1:\", logr.coef_)\n",
    "\n",
    "#ov parametri ne odgovaraju onima koje sam imala pre\n",
    "print(\"w0:\", clf.intercept_)\n",
    "print(\"w1:\", clf.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "G5Pszy_NAtW0",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 484,
     "status": "ok",
     "timestamp": 1733842427359,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "G5Pszy_NAtW0",
    "outputId": "588baab3-2d12-4654-e89a-adf2519fd80d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w0: [-11.68319035]\n",
      "w1: [[5.68708926]]\n"
     ]
    }
   ],
   "source": [
    "logr = linear_model.LogisticRegression(None)\n",
    "clf = logr.fit(X, flags)\n",
    "\n",
    "\n",
    "# underlaying linear model coefficients\n",
    "print(\"w0:\", logr.intercept_)\n",
    "print(\"w1:\", logr.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "54f463d0",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 493
    },
    "executionInfo": {
     "elapsed": 978,
     "status": "ok",
     "timestamp": 1733842571838,
     "user": {
      "displayName": "Marija Delić",
      "userId": "03427129792104028391"
     },
     "user_tz": -60
    },
    "id": "54f463d0",
    "outputId": "a5fad6b1-8c25-4d6f-d7f5-06a740376da5",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x703a6d588530>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 37764 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": 12,
   "id": "2c16aad9",
   "metadata": {
    "id": "2c16aad9",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covariance matrix: \n",
      " [[12.31854544 -6.71118734]\n",
      " [-6.71118734  3.81820749]]\n",
      "Standard errors:  [3.50977855 1.95402341]\n",
      "logitParams: [-11.68319035   5.68708926]\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": 13,
   "id": "2f0ac198",
   "metadata": {
    "id": "2f0ac198",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Logistic regresion modeling with estimated uncertainties')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 46140 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": {
    "id": "6a1a899c-ac8a-44e8-bad6-04262e429913"
   },
   "source": [
    "### Distribution of parameters via BS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "712f33c6-b772-4c90-b2e0-ad9336fc8d83",
   "metadata": {
    "id": "712f33c6-b772-4c90-b2e0-ad9336fc8d83"
   },
   "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": 15,
   "id": "8d8c5f8a-f065-4980-a313-786e0df9b2c2",
   "metadata": {
    "id": "8d8c5f8a-f065-4980-a313-786e0df9b2c2"
   },
   "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": 16,
   "id": "d2ff06ca-e6f3-4442-9bce-29ee929c0071",
   "metadata": {
    "id": "d2ff06ca-e6f3-4442-9bce-29ee929c0071",
    "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",
    "    penalty = 'none' if choice_setup == 0 else None\n",
    "    clf = linear_model.LogisticRegression(penalty = penalty).fit(X[idx], flags[idx])\n",
    "    \n",
    "    return np.insert(clf.coef_, 0, clf.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "69792c6a-18d5-43c1-b1f9-a0744beb1ebb",
   "metadata": {
    "id": "69792c6a-18d5-43c1-b1f9-a0744beb1ebb",
    "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**5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "a97debd2-00cf-4ebc-9b50-da847daf6d63",
   "metadata": {
    "id": "a97debd2-00cf-4ebc-9b50-da847daf6d63",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32360 characters]",
      "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', \n",
    "         verticalalignment='center', transform=ax.figure.transFigure)\n",
    "plt.savefig(os.path.join(figsdir, \"logit_boots\" + figsext), bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4f08a287-fdcb-47cc-a2fd-55ddba066ce8",
   "metadata": {
    "id": "4f08a287-fdcb-47cc-a2fd-55ddba066ce8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size_selectd: 63788\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 43988 characters]",
      "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', \n",
    "         verticalalignment='center', transform=ax.figure.transFigure)\n",
    "plt.plot([p0[0]], [p0[1]], marker=\"o\", markersize=10, markeredgecolor=\"red\", \n",
    "         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(os.path.join(figsdir, \"logit_boots_sel\" + figsext), bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "668f5b02-e351-4dde-9697-af36b49588f2",
   "metadata": {
    "id": "668f5b02-e351-4dde-9697-af36b49588f2"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLE: beta:\n",
      " [-11.68319035   5.68708926]\n",
      "MLE: cov(beta):\n",
      " [[12.31854544 -6.71118734]\n",
      " [-6.71118734  3.81820749]]\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": 21,
   "id": "308b94a4-7a95-47c8-b0c5-4a80bc05b3df",
   "metadata": {
    "id": "308b94a4-7a95-47c8-b0c5-4a80bc05b3df"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample_mean(beta):\n",
      " [-10.94752988   5.03018971]\n",
      "sample_cov(beta):\n",
      " [[ 5.67071502 -2.86912931]\n",
      " [-2.86912931  1.73456244]]\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": 22,
   "id": "cc3d41a7-518f-4804-832a-090cd1501d35",
   "metadata": {
    "id": "cc3d41a7-518f-4804-832a-090cd1501d35",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36152 characters]",
      "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, \n",
    "         horizontalalignment = 'center', \n",
    "         verticalalignment = 'center', \n",
    "         transform = ax.figure.transFigure)\n",
    "plt.plot([p0[0]], [p0[1]], marker=\"o\", \n",
    "         markersize = 10, \n",
    "         markeredgecolor = \"red\",\n",
    "         markerfacecolor = \"red\")\n",
    "\n",
    "plt.xlim(-80,30)\n",
    "plt.ylim(-20,40)\n",
    "plt.savefig(os.path.join(figsdir, \"logit_boots_zoom\" + figsext), bbox_inches='tight')\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4c2f7492-6caa-4f05-8a2f-a5a1eac28d30",
   "metadata": {
    "id": "4c2f7492-6caa-4f05-8a2f-a5a1eac28d30",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact: opt(p)= [-11.68319035   5.68708926]\n",
      "exact: covlogit= [[12.31854544 -6.71118734]\n",
      " [-6.71118734  3.81820749]]\n",
      "mean(p)= [-31.45398028  17.4033474 ]\n",
      "cov(p)= [[ 1728.61003364 -1069.08739753]\n",
      " [-1069.08739753   664.52762393]]\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": 24,
   "id": "9bd0c944-1793-42f5-8a78-8a38e339e8f4",
   "metadata": {
    "id": "9bd0c944-1793-42f5-8a78-8a38e339e8f4",
    "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": 25,
   "id": "5160ac22-29e7-4c06-8b60-1214231fbc43",
   "metadata": {
    "id": "5160ac22-29e7-4c06-8b60-1214231fbc43",
    "tags": []
   },
   "outputs": [],
   "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": {
    "id": "98ad2d95-d89b-4f18-a2f7-c7fc17473135"
   },
   "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": {
    "id": "e065c5d1-7dd5-4098-a095-5a283a8d408c"
   },
   "source": [
    "## Final results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "8305144d",
   "metadata": {
    "id": "8305144d",
    "tags": []
   },
   "outputs": [],
   "source": [
    "def gen_plots(pert, organ, prefix, alpha, title = \"\", \n",
    "              bootstrapping = \"nonparam\", Nsamples = 10**4, \n",
    "              xvert = None):\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",
    "    else:\n",
    "        print(\"This organ is not supported\")\n",
    "        return\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",
    "    penalty = 'none' if choice_setup == 0 else None\n",
    "    clf = linear_model.LogisticRegression(penalty=penalty).fit(X, flags)\n",
    "\n",
    "    # Store fitted parameters in array\n",
    "    logitParams = np.insert(clf.coef_, 0, clf.intercept_)\n",
    "\n",
    "    #print params\n",
    "    print(\"popt:\", logitParams)\n",
    "\n",
    "    # logit function\n",
    "    f = lambda x, p: 1/(1 + np.exp(-p[0] - p[1]*x))\n",
    "\n",
    "    if bootstrapping == \"nonparam\":  # non-param boostrapping \n",
    "        # obtaining parameters via bootstrapping method\n",
    "        rng = np.random.default_rng()\n",
    "\n",
    "        # generate parameters: make a sample by sampling with replacement and fit it\n",
    "        def gen_pars():\n",
    "            idx = rng.integers(len(X), size=len(X))\n",
    "            st = (flags[idx] == 0)\n",
    "            \n",
    "            # a bit constraint sampling with replacement: two groups need to be represented samples\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 = penalty).fit(X[idx], flags[idx])\n",
    "\n",
    "            return np.insert(tmp.coef_, 0, tmp.intercept_)\n",
    "\n",
    "        ps = np.array([gen_pars() for _ in range(Nsamples)])\n",
    "\n",
    "    elif bootstrapping == \"param\" :  # param boostrapping \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",
    "    else:\n",
    "        print(\"This boostrapping  is not supported\")\n",
    "        return\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",
    "\n",
    "    fs = 16\n",
    "    params = {'font.size': fs, 'axes.labelsize': fs,'axes.titlesize':fs, \n",
    "              'legend.fontsize': fs, 'xtick.labelsize': fs, 'ytick.labelsize': fs}\n",
    "    plt.rcParams.update(params)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize = (7, 5))\n",
    "\n",
    "    plt.xlabel(f\"$x = \\\\mathrm{{max}}(SUV(visit,p = {{{pert}}} \\\\%))$\")\n",
    "    plt.ylabel(r\"$P(Y = AE|X = x)$\")\n",
    "   \n",
    "    ax.set_xlabel(f\"$x = \\\\mathrm{{max}}(SUV(visit,p = {{{pert}}} \\\\%))$\")\n",
    "    ax.set_ylabel(f\"$P(Y = AE|X = x)$\")\n",
    "\n",
    "    ax.scatter(X[flags==0], flags[flags==0], label=\"NC\",zorder=2.5)\n",
    "    ax.scatter(X[flags==1], flags[flags==1], label=\"AE\",zorder=2.5)\n",
    "\n",
    "    ax.plot(xs, yopts, label=\"optimal\", color=\"red\",zorder=3)\n",
    "    ax.plot(xs, mys, label=\"median\", color=\"green\",zorder=3)\n",
    "\n",
    "    #plt.fill_between(xs, yerr_min, yerr_max, alpha=0.5, label=\"confidence band\")\n",
    "    ax.fill_between(xs, yerr_min, yerr_max, color='xkcd:beige', alpha=0.7, label=\"confidence interval\")\n",
    "\n",
    "    if title != \"\": ax.set_title(title+\":\" + organ)\n",
    "    if xvert is not None: \n",
    "        ax.axvline(x = xvert, linestyle = \"--\", color = \"k\", alpha = 0.7, label=\"example\")\n",
    "        \n",
    "        I = lambda v: scipy.interpolate.interp1d(xs, v)(xvert)        \n",
    "        display(pd.DataFrame({\"type\": [\"optimal\", \"median\", \"LCL\", \"UCL\"], \n",
    "                              f\"value@{xvert}\": [I(yopts), I(mys), I(yerr_min), I(yerr_max)]}))\n",
    "\n",
    "    ax.text(0.1, 0.93, \"(\" + organ + \")\", fontsize = 16, \n",
    "            horizontalalignment = 'center', verticalalignment = 'center', \n",
    "            transform = fig.transFigure)\n",
    "\n",
    "    ax.legend(loc = 'center right', fancybox=True, fontsize='10', framealpha=0.5)\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    \n",
    "    basename =  f\"{prefix}_{organ}{figsext}\"\n",
    "    plt.savefig(os.path.join(figsdir, basename))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb7dda9a-2deb-44ac-823d-66773165be59",
   "metadata": {
    "id": "eb7dda9a-2deb-44ac-823d-66773165be59"
   },
   "source": [
    "### Parametric boostrapping : MVN approx of params (! primary choice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "593f4397-ea4b-44b4-8249-0502f8681ce7",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "3b3a0f89-bad4-4a39-81b9-1b9ac8ca3f75",
   "metadata": {
    "id": "3b3a0f89-bad4-4a39-81b9-1b9ac8ca3f75",
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-11.68319035   5.68708926]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>type</th>\n",
       "      <th>value@2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>optimal</td>\n",
       "      <td>0.4233901220524278</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>median</td>\n",
       "      <td>0.42470951308011934</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LCL</td>\n",
       "      <td>0.11738712465726915</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>UCL</td>\n",
       "      <td>0.8032612583732215</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      type              value@2\n",
       "0  optimal   0.4233901220524278\n",
       "1   median  0.42470951308011934\n",
       "2      LCL  0.11738712465726915\n",
       "3      UCL   0.8032612583732215"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 61832 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-3.47422933  0.33708866]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 59240 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-4.3866421   1.16301187]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 61704 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12.8 s, sys: 2.21 s, total: 15 s\n",
      "Wall time: 12.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(95, \"lung\", \"logit_mvn\", alpha, bootstrapping=\"param\", xvert=2)\n",
    "gen_plots(95, \"bowel\", \"logit_mvn\", alpha, bootstrapping=\"param\")\n",
    "gen_plots(75, \"thyroid\", \"logit_mvn\", alpha, bootstrapping=\"param\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1e9c470-9a7d-4490-9912-be9d24a31bab",
   "metadata": {
    "id": "b1e9c470-9a7d-4490-9912-be9d24a31bab"
   },
   "source": [
    "### Non-parametric bootstrapping: sampling with replacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "1ceb9959-7d52-4821-a07f-22cbaafc019a",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "2c461077-18cd-4200-88fb-f10976449eb6",
   "metadata": {
    "id": "2c461077-18cd-4200-88fb-f10976449eb6"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-11.68319035   5.68708926]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 59856 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-3.47422933  0.33708866]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 62400 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-4.3866421   1.16301187]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_199323/3990547804.py:34: RuntimeWarning: overflow encountered in exp\n",
      "  f = lambda x, p: 1/(1 + np.exp(-p[0] - p[1]*x))\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 67460 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 25s, sys: 1.95 s, total: 1min 27s\n",
      "Wall time: 1min 25s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(95, \"lung\", \"logit_bs\", alpha, bootstrapping = \"nonparam\")\n",
    "gen_plots(95, \"bowel\", \"logit_bs\", alpha, bootstrapping = \"nonparam\")\n",
    "gen_plots(75, \"thyroid\", \"logit_bs\", alpha, bootstrapping = \"nonparam\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20503b2a-51c4-4ad1-8cff-d222d734e198",
   "metadata": {
    "id": "20503b2a-51c4-4ad1-8cff-d222d734e198"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
