{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "078900c6-d71d-44f0-9e76-cc5ce66f0ca9",
   "metadata": {},
   "source": [
    "# Logistic regression: fit\n",
    "\n",
    "Exploring different fitting techniques and using MLE asymptotic approx of parameter distribution to estimate model CI intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e18cec5e",
   "metadata": {},
   "source": [
    "## Common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3b351af7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "from sklearn import linear_model\n",
    "\n",
    "# our libs\n",
    "import data_utils\n",
    "import logit_utils"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0de7f02-942b-42ff-874c-3c0e2ac1e0dc",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3373f87f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_path = \"./results\"\n",
    "if not os.path.exists(results_path): os.makedirs(results_path)\n",
    "\n",
    "data_path = \"../../data/\"\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": 3,
   "id": "811d990d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Loading all data\n",
    "\"\"\"\n",
    "suv_dict = scipy.io.loadmat(suv_filename)\n",
    "flags_dict = scipy.io.loadmat(flags_filename)\n",
    "normal_range_dict = scipy.io.loadmat(normal_range_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4a292c31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get data\n",
    "perc = 95\n",
    "organ = \"lung\"\n",
    "\n",
    "x0, y = data_utils.get_data(organ, perc, suv_dict, flags_dict)\n",
    "\n",
    "labs = [\"plain\", \"log\"]\n",
    "xs = [x0, np.log(x0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4",
   "metadata": {},
   "source": [
    "## Fit\n",
    "\n",
    "Ref: 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": 5,
   "id": "239ae789",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parss = array([[-11.68319035,   5.68708926],\n",
      "       [ -7.52531005,  10.73177719]])\n"
     ]
    }
   ],
   "source": [
    "# fits\n",
    "logr = linear_model.LogisticRegression(penalty = None)\n",
    "parss = np.array([logit_utils.logit_poly_fit(logr, x, y, degree = 1) for x in xs])\n",
    "print(f\"{parss = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9cf6fec1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cov_parss = array([[[12.31854544, -6.71118734],\n",
      "        [-6.71118734,  3.81820749]],\n",
      "\n",
      "       [[ 5.04803783, -7.92922372],\n",
      "        [-7.92922372, 13.99363052]]])\n"
     ]
    }
   ],
   "source": [
    "# asymptotic covariance matrix of parameters\n",
    "cov_parss = np.array([logit_utils.logit_poly_cov(x, pars) for x, pars in zip(xs, parss)])\n",
    "print(f\"{cov_parss = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2359d3c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining two-sided confidence intervals\n",
    "alpha = 0.05   # significance level\n",
    "probs = [alpha/2, 1 - alpha/2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4ba4d5d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pars_CIs = array([[[-18.5622299 ,   1.85727375],\n",
      "        [ -4.80415081,   9.51690477]],\n",
      "\n",
      "       [[-11.92892555,   3.3999319 ],\n",
      "        [ -3.12169455,  18.06362248]]])\n"
     ]
    }
   ],
   "source": [
    "# CI of params (assuming asymptotic distr of parameters)\n",
    "pars_CIs = np.array([logit_utils.logit_poly_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
    "print(f\"{pars_CIs = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "96acc3ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "('scale', 'coef')",
         "rawType": "object",
         "type": "unknown"
        },
        {
         "name": "value",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "LCL",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "UCL",
         "rawType": "float64",
         "type": "float"
        }
       ],
       "ref": "d7e2eb90-ee4f-4a41-9d9a-059fe24771ed",
       "rows": [
        [
         "('plain', 'b0')",
         "-11.683190353932442",
         "-18.562229897321515",
         "-4.804150810543373"
        ],
        [
         "('plain', 'b1')",
         "5.687089258005262",
         "1.8572737464039784",
         "9.516904769606544"
        ],
        [
         "('log', 'b0')",
         "-7.5253100481979",
         "-11.928925546652636",
         "-3.1216945497431663"
        ],
        [
         "('log', 'b1')",
         "10.731777187587987",
         "3.399931895837665",
         "18.063622479338306"
        ]
       ],
       "shape": {
        "columns": 3,
        "rows": 4
       }
      },
      "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></th>\n",
       "      <th>value</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>scale</th>\n",
       "      <th>coef</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">plain</th>\n",
       "      <th>b0</th>\n",
       "      <td>-11.683190</td>\n",
       "      <td>-18.562230</td>\n",
       "      <td>-4.804151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1</th>\n",
       "      <td>5.687089</td>\n",
       "      <td>1.857274</td>\n",
       "      <td>9.516905</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">log</th>\n",
       "      <th>b0</th>\n",
       "      <td>-7.525310</td>\n",
       "      <td>-11.928926</td>\n",
       "      <td>-3.121695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1</th>\n",
       "      <td>10.731777</td>\n",
       "      <td>3.399932</td>\n",
       "      <td>18.063622</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                value        LCL        UCL\n",
       "scale coef                                 \n",
       "plain b0   -11.683190 -18.562230  -4.804151\n",
       "      b1     5.687089   1.857274   9.516905\n",
       "log   b0    -7.525310 -11.928926  -3.121695\n",
       "      b1    10.731777   3.399932  18.063622"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# present table of results\n",
    "d = parss.shape[-1]\n",
    "index = pd.MultiIndex.from_tuples([(lab, f\"b{i}\") for lab in labs for i in range(d)], names = ['scale', 'coef'])\n",
    "pd.DataFrame({\"value\": parss.flatten(), \"LCL\": pars_CIs[:,0].flatten(), \"UCL\": pars_CIs[:,1].flatten()}, index = index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6ef6f46b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "scale",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "LLF",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "AIC",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "BIC",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "chi2",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "p-value(chi2)",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "A",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "n",
         "rawType": "int64",
         "type": "integer"
        },
        {
         "name": "k",
         "rawType": "int64",
         "type": "integer"
        },
        {
         "name": "dof",
         "rawType": "int64",
         "type": "integer"
        }
       ],
       "ref": "41b39871-5b55-41d6-a947-656a9aa52de7",
       "rows": [
        [
         "plain",
         "-7.394195827272766",
         "18.78839165454553",
         "22.90927767563837",
         "23.105937357608315",
         "0.9999705652971257",
         "0.9655172413793104",
         "58",
         "2",
         "56"
        ],
        [
         "log",
         "-6.795378655860808",
         "17.590757311721617",
         "21.711643332814454",
         "19.34046691560593",
         "0.9999987944563601",
         "0.9655172413793104",
         "58",
         "2",
         "56"
        ]
       ],
       "shape": {
        "columns": 9,
        "rows": 2
       }
      },
      "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>LLF</th>\n",
       "      <th>AIC</th>\n",
       "      <th>BIC</th>\n",
       "      <th>chi2</th>\n",
       "      <th>p-value(chi2)</th>\n",
       "      <th>A</th>\n",
       "      <th>n</th>\n",
       "      <th>k</th>\n",
       "      <th>dof</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>scale</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>plain</th>\n",
       "      <td>-7.394196</td>\n",
       "      <td>18.788392</td>\n",
       "      <td>22.909278</td>\n",
       "      <td>23.105937</td>\n",
       "      <td>0.999971</td>\n",
       "      <td>0.965517</td>\n",
       "      <td>58</td>\n",
       "      <td>2</td>\n",
       "      <td>56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log</th>\n",
       "      <td>-6.795379</td>\n",
       "      <td>17.590757</td>\n",
       "      <td>21.711643</td>\n",
       "      <td>19.340467</td>\n",
       "      <td>0.999999</td>\n",
       "      <td>0.965517</td>\n",
       "      <td>58</td>\n",
       "      <td>2</td>\n",
       "      <td>56</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            LLF        AIC        BIC       chi2  p-value(chi2)         A   n  \\\n",
       "scale                                                                           \n",
       "plain -7.394196  18.788392  22.909278  23.105937       0.999971  0.965517  58   \n",
       "log   -6.795379  17.590757  21.711643  19.340467       0.999999  0.965517  58   \n",
       "\n",
       "       k  dof  \n",
       "scale          \n",
       "plain  2   56  \n",
       "log    2   56  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# goodness of fit measures\n",
    "index = pd.Index(labs, name = 'scale')\n",
    "pd.DataFrame([logit_utils.logit_poly_goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], index = index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "54f463d0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 94996 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting results\n",
    "fig, axs = plt.subplots(ncols = len(labs), figsize = (6*len(labs), 5))\n",
    "\n",
    "for ax, lab, x, pars, cov_pars in zip(axs, labs, xs, parss, cov_parss):\n",
    "\n",
    "    ax.set_title(f\"cond. prob. for AE: scale {lab}\")\n",
    "\n",
    "    if lab == \"plain\":\n",
    "        xlab = r\"$x = max_{visit} SUV(visit, p)$\"\n",
    "    else:\n",
    "        xlab = r\"$x = log(max_{visit} SUV(visit, p))$\"\n",
    "    \n",
    "    ax.set_xlabel(xlab)\n",
    "    ax.set_ylabel(\"P(AE|X = x)\")\n",
    "\n",
    "    # plotting data points\n",
    "    for i, lab in enumerate([\"NC\", \"AE\"]): ax.scatter(x[y==i], y[y==i], label = lab)  \n",
    "\n",
    "    # plotting fitted model \n",
    "    xp = np.linspace(min(x), max(x), 100)\n",
    "    yp = logit_utils.logit_poly_model(xp, pars)\n",
    "\n",
    "    ax.plot(xp, yp, label = \"logistic reg\")\n",
    "\n",
    "    for lab, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
    "        \n",
    "        # define quantile model\n",
    "        fname = f\"logit_utils.logit_poly_model_quantiles_{lab}\"\n",
    "        \n",
    "        # get result of model quantiles at probs\n",
    "        # pars and cov_pars are obtained via MLE method\n",
    "        res = eval(fname)(xp, probs, pars, cov_pars)\n",
    "\n",
    "        # make plot of quantiles\n",
    "        ax.fill_between(xp, *res, color = c, alpha = 0.5, label = f\"CI:{lab}\")\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "fig.legend(handles, labels, bbox_to_anchor=(1.02, 0.5), loc='center right')\n",
    "\n",
    "plt.savefig(os.path.join(results_path, \"logit_fit.pdf\"))\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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
}
