{
"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": 2,
"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",
"# 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": 3,
"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": 4,
"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": 5,
"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": "code",
"execution_count": 6,
"id": "130de003",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "7d1554b0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0 1 2 3 4]\n",
" [5 6 7 8 9]]\n",
"[[[ 0 1 2]\n",
" [ 3 4 5]\n",
" [ 6 7 8]\n",
" [ 9 10 11]\n",
" [12 13 14]]\n",
"\n",
" [[15 16 17]\n",
" [18 19 20]\n",
" [21 22 23]\n",
" [24 25 26]\n",
" [27 28 29]]]\n"
]
}
],
"source": [
"a=np.arange(10).reshape(2,5)\n",
"\n",
"b=np.arange(10*3).reshape(2,5,3)\n",
"\n",
"print(a)\n",
"print(b)"
]
},
{
"cell_type": "markdown",
"id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4",
"metadata": {},
"source": [
"## Fit linear model\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": 8,
"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",
"parss = np.array([logit_utils.logit_poly_fit(x, y, degree = 1) for x in xs])\n",
"print(f\"{parss = }\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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": 9,
"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": "193de933-00a7-456e-bba7-b455167036d8",
"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": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" value | \n",
" LCL | \n",
" UCL | \n",
"
\n",
" \n",
" | scale | \n",
" coef | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | plain | \n",
" b0 | \n",
" -11.683190 | \n",
" -18.562230 | \n",
" -4.804151 | \n",
"
\n",
" \n",
" | b1 | \n",
" 5.687089 | \n",
" 1.857274 | \n",
" 9.516905 | \n",
"
\n",
" \n",
" | log | \n",
" b0 | \n",
" -7.525310 | \n",
" -11.928926 | \n",
" -3.121695 | \n",
"
\n",
" \n",
" | b1 | \n",
" 10.731777 | \n",
" 3.399932 | \n",
" 18.063622 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" LLF | \n",
" AIC | \n",
" BIC | \n",
" A | \n",
" chi2 | \n",
" p-value(chi2) | \n",
" n | \n",
" k | \n",
" dof | \n",
"
\n",
" \n",
" | scale | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | plain | \n",
" -7.394196 | \n",
" 18.788392 | \n",
" 22.909278 | \n",
" 0.965517 | \n",
" 23.105937 | \n",
" 0.999971 | \n",
" 58 | \n",
" 2 | \n",
" 56 | \n",
"
\n",
" \n",
" | log | \n",
" -6.795379 | \n",
" 17.590757 | \n",
" 21.711643 | \n",
" 0.965517 | \n",
" 19.340467 | \n",
" 0.999999 | \n",
" 58 | \n",
" 2 | \n",
" 56 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"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"
]
},
{
"cell_type": "markdown",
"id": "50dfb646",
"metadata": {},
"source": [
"## Fit cubic model\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": 2,
"id": "bcb24955",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "855f597f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([10, 13, 16, 19, 22])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.array([1,2])@np.arange(10).reshape(2,5)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "896ca504",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"parss = array([[-197.73509939, 248.907034 , -98.78754738, 12.20769479],\n",
" [ -75.19604143, 286.98532686, -339.07257923, 124.33155901]])\n"
]
}
],
"source": [
"# fits\n",
"parss = np.array([logit_utils.logit_poly_fit(x, y, degree = 3) for x in xs])\n",
"print(f\"{parss = }\")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "45404e51",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"cov_parss = array([[[ 11032.93517532, -14145.10251152, 5717.53233434,\n",
" -718.21924319],\n",
" [-14145.10251152, 18209.22228607, -7393.24313207,\n",
" 933.58003493],\n",
" [ 5717.53233434, -7393.24313207, 3017.99869597,\n",
" -383.74314053],\n",
" [ -718.21924319, 933.58003493, -383.74314053,\n",
" 49.27459665]],\n",
"\n",
" [[ 1992.00030798, -7682.85412745, 9175.28509088,\n",
" -3379.47618792],\n",
" [ -7682.85412745, 29867.99665855, -35928.29336062,\n",
" 13325.79628235],\n",
" [ 9175.28509088, -35928.29336062, 43549.830978 ,\n",
" -16289.96884834],\n",
" [ -3379.47618792, 13325.79628235, -16289.96884834,\n",
" 6159.65134352]]])\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": 14,
"id": "4b2574bf",
"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": 15,
"id": "f3f7a0cc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pars_CIs = array([[[-403.60536472, -15.57354671, -206.46074605, -1.55044222],\n",
" [ 8.13516594, 513.3876147 , 8.88565129, 25.9658318 ]],\n",
"\n",
" [[-162.67282209, -51.74270449, -748.08954997, -29.49316913],\n",
" [ 12.28073922, 625.71335821, 69.94439151, 278.15628715]]])\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": 16,
"id": "ad95bf42",
"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": "3591ba44-5738-4dfa-9dc9-648bcbda1322",
"rows": [
[
"('plain', 'b0')",
"-197.73509939400154",
"-403.60536472485205",
"8.13516593684895"
],
[
"('plain', 'b1')",
"248.90703399728125",
"-15.57354670535642",
"513.3876146999189"
],
[
"('plain', 'b2')",
"-98.78754738088045",
"-206.460746049195",
"8.885651287434086"
],
[
"('plain', 'b3')",
"12.20769479113704",
"-1.550442217446367",
"25.965831799720444"
],
[
"('log', 'b0')",
"-75.19604143281751",
"-162.67282208873257",
"12.280739223097527"
],
[
"('log', 'b1')",
"286.98532685781333",
"-51.74270448955622",
"625.7133582051829"
],
[
"('log', 'b2')",
"-339.07257923019125",
"-748.0895499663287",
"69.94439150594616"
],
[
"('log', 'b3')",
"124.33155901016768",
"-29.4931691252947",
"278.15628714563"
]
],
"shape": {
"columns": 3,
"rows": 8
}
},
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" value | \n",
" LCL | \n",
" UCL | \n",
"
\n",
" \n",
" | scale | \n",
" coef | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | plain | \n",
" b0 | \n",
" -197.735099 | \n",
" -403.605365 | \n",
" 8.135166 | \n",
"
\n",
" \n",
" | b1 | \n",
" 248.907034 | \n",
" -15.573547 | \n",
" 513.387615 | \n",
"
\n",
" \n",
" | b2 | \n",
" -98.787547 | \n",
" -206.460746 | \n",
" 8.885651 | \n",
"
\n",
" \n",
" | b3 | \n",
" 12.207695 | \n",
" -1.550442 | \n",
" 25.965832 | \n",
"
\n",
" \n",
" | log | \n",
" b0 | \n",
" -75.196041 | \n",
" -162.672822 | \n",
" 12.280739 | \n",
"
\n",
" \n",
" | b1 | \n",
" 286.985327 | \n",
" -51.742704 | \n",
" 625.713358 | \n",
"
\n",
" \n",
" | b2 | \n",
" -339.072579 | \n",
" -748.089550 | \n",
" 69.944392 | \n",
"
\n",
" \n",
" | b3 | \n",
" 124.331559 | \n",
" -29.493169 | \n",
" 278.156287 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" LLF | \n",
" AIC | \n",
" BIC | \n",
" A | \n",
" chi2 | \n",
" p-value(chi2) | \n",
" n | \n",
" k | \n",
" dof | \n",
"
\n",
" \n",
" | scale | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | plain | \n",
" -3.188317 | \n",
" 14.376634 | \n",
" 22.618406 | \n",
" 0.965517 | \n",
" 6.024924 | \n",
" 1.0 | \n",
" 58 | \n",
" 4 | \n",
" 54 | \n",
"
\n",
" \n",
" | log | \n",
" -3.491426 | \n",
" 14.982853 | \n",
" 23.224625 | \n",
" 0.965517 | \n",
" 6.736954 | \n",
" 1.0 | \n",
" 58 | \n",
" 4 | \n",
" 54 | \n",
"
\n",
" \n",
"
\n",
"