{
 "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",
    "    * stratified non-parametric bootstrapping\n",
    "    * parametric bootstrapping of parameters around highest peak obtained via non-parametric bootstrapping or global optimimum"
   ]
  },
  {
   "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",
    "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_sb\"\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 0x7b891f5e9880>"
      ]
     },
     "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, 45524 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, 30200 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: 63533\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 44136 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.94139483   5.02655264]\n",
      "sample_cov(beta):\n",
      " [[ 5.55949878 -2.84410312]\n",
      " [-2.84410312  1.73167577]]\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, 35844 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.69607061  17.54685413]\n",
      "cov(p)= [[ 1748.98545671 -1081.72138611]\n",
      " [-1081.72138611   672.4002813 ]]\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": [
    "\n",
    "\n",
    "def gen_plots(pert, organ, prefix, alpha, title = \"\", bootstrapping = \"nonparam\", Nsamples = 10**3):\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",
    "    Y = 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, Y)\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",
    "        Xs = [X[Y==i] for i in range(2)]\n",
    "        Ns = [len(Xs[i]) for i in range(2)]\n",
    "\n",
    "        ps = []\n",
    "        for _ in range(Nsamples):\n",
    "            # to the groups stratified sampling\n",
    "            \n",
    "            # stratified sampling with replacement\n",
    "            Xb = [rng.choice(Xs[i], Ns[i]) for i in range(2)]\n",
    "            Yb = [np.full(Ns[i], i) for i in range(2)]\n",
    "            \n",
    "            Xb = np.concatenate(Xb)\n",
    "            Yb = np.concatenate(Yb)\n",
    "\n",
    "            tmp = linear_model.LogisticRegression(penalty = penalty).fit(Xb, Yb)\n",
    "            pfit = np.insert(tmp.coef_, 0, tmp.intercept_)\n",
    "        \n",
    "            ps.append(pfit)   \n",
    "        ps = np.array(ps)\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",
    "    plt.clf()\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 = plt.figure(figsize = (7, 5))\n",
    "        \n",
    "    plt.xlabel(f\"$x = \\\\mathrm{{max}}(SUV(visit,p = {{{pert}}} \\\\%))$\")\n",
    "    plt.ylabel(f\"$P(Y = AE|X = x)$\")\n",
    "\n",
    "    plt.scatter(X[flags==0], flags[flags==0], label=\"NC\",zorder=2.5)\n",
    "    plt.scatter(X[flags==1], flags[flags==1], label=\"AE\",zorder=2.5)\n",
    "\n",
    "    plt.plot(xs, yopts, label=\"optimal\", color=\"red\",zorder=3)\n",
    "    plt.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",
    "    plt.fill_between(xs, yerr_min, yerr_max, color='xkcd:beige', alpha=0.5, label=\"prediction interval\")\n",
    "\n",
    "    if title != \"\": plt.title(title+\":\" + organ)\n",
    "   \n",
    "    plt.text(0.1, 0.93, \"(\" + organ + \")\", fontsize=16, horizontalalignment='center', \n",
    "             verticalalignment='center', transform=fig.transFigure)\n",
    "\n",
    "    plt.legend(loc = 'center right', fancybox=True, fontsize='10', framealpha=0.5)\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    basename =  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",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "popt: [-11.68319035   5.68708926]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 58352 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": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 61356 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": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 63036 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.27 s, sys: 2.35 s, total: 5.62 s\n",
      "Wall time: 2.89 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gen_plots(95, \"lung\", \"logit_mvn\", alpha, bootstrapping=\"param\")\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": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 58496 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": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 62228 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_23628/878945397.py:32: RuntimeWarning: overflow encountered in exp\n",
      "  f = lambda x, p: 1/(1 + np.exp(-p[0] - p[1]*x))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 66344 characters]",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.3 s, sys: 1.86 s, total: 12.1 s\n",
      "Wall time: 9.93 s\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
}
