{
"cells": [
{
"cell_type": "markdown",
"id": "be915d56-9a8a-422b-9da4-d7ca37dd0b85",
"metadata": {
"id": "be915d56-9a8a-422b-9da4-d7ca37dd0b85"
},
"source": [
"# Bayesian modelling\n",
"\n",
"Two separate-group discussion: \n",
"\n",
" 1. fitting each groups separately \n",
" 2. obtained distributions of each group are then used in conditional probability to encounter measurement in the second group."
]
},
{
"cell_type": "markdown",
"id": "6fa3ccb1-05c0-4ab9-a74a-a4e20d1e38f8",
"metadata": {
"id": "6fa3ccb1-05c0-4ab9-a74a-a4e20d1e38f8"
},
"source": [
"## Data & common"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ee1ce969-7f96-4ca9-882e-866805bd316c",
"metadata": {},
"outputs": [],
"source": [
"choice_setup = 1 # 0 for MD, 1 for MH"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "uq21BzGM7YHZ",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"executionInfo": {
"elapsed": 26515,
"status": "ok",
"timestamp": 1715067825951,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "uq21BzGM7YHZ",
"outputId": "fed613c4-2be9-49ba-bdc5-0a09399dc2fd"
},
"outputs": [],
"source": [
"if choice_setup == 0:\n",
" from google.colab import drive\n",
" drive.mount('/content/drive')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0de55416",
"metadata": {
"executionInfo": {
"elapsed": 1711,
"status": "ok",
"timestamp": 1715067832083,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "0de55416",
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import pandas as pd\n",
"\n",
"import scipy\n",
"from scipy import io\n",
"from scipy import stats\n",
"from scipy.optimize import minimize\n",
"\n",
"import os"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "71b466ba",
"metadata": {
"executionInfo": {
"elapsed": 305,
"status": "ok",
"timestamp": 1715067835069,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "71b466ba",
"tags": []
},
"outputs": [],
"source": [
"wd = os.getcwd() # working directory\n",
"\n",
"topic = \"bayes\"\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}/\" # figure directory\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",
"\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": "c72c613c",
"metadata": {
"executionInfo": {
"elapsed": 2756,
"status": "ok",
"timestamp": 1715067840744,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "c72c613c",
"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": "8c4a5c84-c72a-4b48-a92c-7e25f5b752fc",
"metadata": {
"id": "8c4a5c84-c72a-4b48-a92c-7e25f5b752fc"
},
"source": [
"## Results"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "eeeaa8fe",
"metadata": {
"executionInfo": {
"elapsed": 11,
"status": "ok",
"timestamp": 1715067844510,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "eeeaa8fe",
"tags": []
},
"outputs": [],
"source": [
"def myfitter(data, fun, cov = True):\n",
"\n",
" nllf = lambda p: -fun(data, *p).sum()\n",
"\n",
" p0 = (1, 0.6)\n",
" bounds= ((0, 10), (0.1, 4))\n",
" res = minimize(nllf, p0, bounds = bounds, method=\"L-BFGS-B\")\n",
"\n",
" if not cov: return res.x\n",
"\n",
" return res.x, res.hess_inv.todense()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "8262d64c",
"metadata": {
"executionInfo": {
"elapsed": 528,
"status": "ok",
"timestamp": 1715067865707,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "8262d64c",
"tags": []
},
"outputs": [],
"source": [
"def gen_plots(pert, organ, dist, 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",
" if dist == \"trunc_Gauss\":\n",
" # create a custom function (\"truncated\" Gaussian distribution)\n",
" fun = lambda x, mu, sigma: -(x-mu)**2/(2*sigma**2) + np.log(2/np.pi)/2 - np.log(sigma*(1 + scipy.special.erf(mu/(np.sqrt(2)*sigma))))\n",
" g = lambda x, mu, sigma: np.exp(-(x-mu)**2/(2*sigma**2))*np.sqrt(2/np.pi)/(sigma*(1 + scipy.special.erf(mu/(np.sqrt(2)*sigma))))\n",
" elif dist == \"Gauss\":\n",
" fun = lambda x, mu, sigma: -(x-mu)**2/(2*sigma**2) - np.log(sigma*np.sqrt(2*np.pi))\n",
" g = lambda x, mu, sigma: np.exp(-(x-mu)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))\n",
" elif dist == \"log_normal\":\n",
" fun = lambda x, mu, sigma: -(np.log(x) - mu)**2/(2*sigma**2) - np.log(x*sigma*np.sqrt(2*np.pi))\n",
" g = lambda x, mu, sigma: np.exp(-(np.log(x)-mu)**2/(2*sigma**2))/(x*sigma*np.sqrt(2*np.pi))\n",
" else:\n",
" print(\"This distribution 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",
" # extracting X associated to NC and AE\n",
" X_NC = X[flags == 0]\n",
" X_AE = X[flags == 1]\n",
" \n",
" #and flags associated to NC and AE\n",
" flags_NC = flags[flags==0]\n",
" flags_AE = flags[flags==1]\n",
"\n",
" popt_NC, pcov_NC = myfitter(X_NC, fun)\n",
" popt_AE, pcov_AE = myfitter(X_AE, fun)\n",
"\n",
" # print optimization results\n",
" print(\"popt_NC:\", popt_NC)\n",
" print(\"popt_AE:\", popt_AE)\n",
"\n",
" print(\"pcov_NC:\", pcov_NC)\n",
" print(\"pcov_AE:\", pcov_AE)\n",
"\n",
" popt = np.concatenate((popt_NC, popt_AE))\n",
"\n",
" if bootstrapping == \"nonparam\":\n",
" # obtaining parameters via bootstrapping method\n",
" rng = np.random.default_rng()\n",
"\n",
" # helper function to sample with replacement of same size sample\n",
" R = lambda X: X[rng.integers(len(X), size=len(X))]\n",
"\n",
" ps = np.array([np.concatenate((myfitter(R(X_NC), fun, False), \n",
" myfitter(R(X_AE), fun, False))) \n",
" for _ in range(Nsamples)])\n",
" elif bootstrapping == \"param\":\n",
" # approximationg parameters with multivatiate distribution\n",
" Z = np.zeros((len(popt_NC), len(popt_AE)))\n",
" CovMatrix = np.bmat([[pcov_NC, Z], [Z, pcov_AE]])\n",
"\n",
" #generate parameters and points in x axis\n",
" ps = np.random.multivariate_normal(popt, CovMatrix, Nsamples)\n",
" else:\n",
" print(\"This boostrapping is not supported\")\n",
" return\n",
"\n",
" \n",
" # correct for outliers\n",
" mask = (ps[:,1]>0) & (ps[:,3]>0)\n",
" ps = ps[mask]\n",
"\n",
" print(\"faulty:\", np.count_nonzero(~mask)/len(mask))\n",
"\n",
" # probability for NC and AE\n",
" P_NC = np.count_nonzero(flags == 0)/len(flags) # Prob(Y = NC)\n",
" P_AE = np.count_nonzero(flags == 1)/len(flags) # Prob(Y = AE)\n",
"\n",
" #firstly for optimal parameter values\n",
" f_XatNC_opt = lambda x: g(x, popt_NC[0], popt_NC[1]) # f(x | NC)\n",
" f_XatAE_opt = lambda x: g(x, popt_AE[0], popt_AE[1]) # f(x | AE)\n",
"\n",
" # probability density for X \n",
" f_X_opt = lambda x: P_NC*f_XatNC_opt(x) + P_AE*f_XatAE_opt(x)\n",
"\n",
" # conditional probability for AE at X\n",
" P_AEatX_opt = lambda x: P_AE*f_XatAE_opt(x)/f_X_opt(x)\n",
"\n",
" # conditional probability for AE at X for generated parameters\n",
" P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
"\n",
" #generate parameters and points in x axis\n",
" xs = np.linspace(1e-3, X.max()+1, 200)\n",
"\n",
" # generate values in y axis\n",
" ys = np.array([np.fromiter(map(lambda p: P_AEatX(x, p), ps),float) for x in xs])\n",
" yopts = P_AEatX_opt(xs)\n",
"\n",
" # doing statistics on values\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",
" 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",
" 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}_{dist}_{organ}{figsext}\"\n",
" plt.savefig(os.path.join(figsdir,basename))\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"id": "a974ab00-e992-4de8-8a9f-431264fe376a",
"metadata": {
"id": "a974ab00-e992-4de8-8a9f-431264fe376a"
},
"source": [
"### Non-parametric boostrapping: sampling with replacement (primary choice)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "00ea6e5e-5e09-45f1-a8e3-0499261ac781",
"metadata": {},
"outputs": [],
"source": [
"alpha = 0.05"
]
},
{
"cell_type": "markdown",
"id": "e79e7dc7-c3b3-46bc-a11b-fc44c3b2f305",
"metadata": {},
"source": [
"#### Lung"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "BDoKoZ0dMekG",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 1000
},
"executionInfo": {
"elapsed": 75311,
"status": "ok",
"timestamp": 1715067952192,
"user": {
"displayName": "Marija Delić",
"userId": "03427129792104028391"
},
"user_tz": -120
},
"id": "BDoKoZ0dMekG",
"outputId": "4498b0d8-e711-414b-9f12-2ea62281d491"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"popt_NC: [1.30761056 0.25845619]\n",
"popt_AE: [2.45119943 0.8708103 ]\n",
"pcov_NC: [[1.18087018e-03 1.19806173e-06]\n",
" [1.19806173e-06 6.30252819e-04]]\n",
"pcov_AE: [[0.16908412 0.03501477]\n",
" [0.03501477 0.07755902]]\n",
"faulty: 0.0\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" type | \n",
" value@2 | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" optimal | \n",
" 0.4699474053976193 | \n",
"
\n",
" \n",
" | 1 | \n",
" median | \n",
" 0.6400362738480477 | \n",
"
\n",
" \n",
" | 2 | \n",
" LCL | \n",
" 0.12335997081396466 | \n",
"
\n",
" \n",
" | 3 | \n",
" UCL | \n",
" 0.9945968429081302 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"popt_NC: [0.25012563 0.18862822]\n",
"popt_AE: [0.84495162 0.3054373 ]\n",
"pcov_NC: [[ 6.83750245e-04 -8.23179537e-06]\n",
" [-8.23179537e-06 3.41031588e-04]]\n",
"pcov_AE: [[0.02184912 0.00118759]\n",
" [0.00118759 0.00977597]]\n",
"faulty: 0.0\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in exp\n",
" P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n",
"/tmp/ipykernel_204016/1182414770.py:100: RuntimeWarning: overflow encountered in scalar multiply\n",
" P_AEatX = lambda x, p: 1/(1 + P_NC/P_AE*np.exp(fun(x, p[0], p[1])-fun(x, p[2], p[3])))\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" type | \n",
" value@2 | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" optimal | \n",
" 0.4482002564354698 | \n",
"
\n",
" \n",
" | 1 | \n",
" median | \n",
" 0.5427501460938366 | \n",
"
\n",
" \n",
" | 2 | \n",
" LCL | \n",
" 0.1843888655804964 | \n",
"
\n",
" \n",
" | 3 | \n",
" UCL | \n",
" 0.9213506594067791 | \n",
"
\n",
" \n",
"
\n",
"