qra/paper_MD/risk_modeling.tex

\documentclass[11pt]{iopart}

% a hack for using amsmath & iopart
\expandafter\let\csname equation*\endcsname\relax
\expandafter\let\csname endequation*\endcsname\relax
\usepackage{amsmath}

\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{subfig}
\usepackage{hyperref}
\usepackage{siunitx}
\usepackage{array}  % For creating new column types


\captionsetup[subfigure]{
  position=top,
  captionskip=-1em,
  singlelinecheck=false,
}
\graphicspath{{./figs/}}

\bibliographystyle{iopart-num}

\newcommand{\nc}{\mathrm{NC}}
\renewcommand{\ae}{\mathrm{AE}}
\newcommand{\ave}[1]{\left\langle #1 \right\rangle}

% Author's guide:
%  https://publishingsupport.iopscience.iop.org/journals/physics-in-medicine-biology/
%

\begin{document}

\title[Quantifying Risk Assessment in Immune Checkpoint Inhibitors]{Quantifying Risk Assessment in Immune Checkpoint Inhibitors Treatment of Malignant Melanoma}
%[Modelling the probability of toxicity with discussion of uncertainties]
%{Modelling the probability of toxicity with the discussion of uncertainties in the therapy treatment of malignant melanoma}

\author{Marija Delić$^{2}$, Martin Horvat$^{1}$, Katja Strašek$^1$, Daniel Huff$^{3}$, Robert Jeraj$^{1,3}$}

\address{%
$^1$Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana, Jadranska cesta 19, 1000 Ljubljana, SI\\
$^2$ Faculty of Technical Sciences, University of Novi Sad, Trg Dositeja Obradovića 6, 21102 Novi Sad, SR \\
$^3$ Department of Medical Physics, School of Medicine and Public Health, University of Wisconsin-Madison, Madison, WI, USA}
\ead{marijadelic@uns.ac.rs}

\begin{abstract}
%
Immune Checkpoint Inhibitors (ICI) have significantly improved cancer treatment outcomes but are associated with a substantial risk of immune-related adverse events (irAEs). Accurate quantification of this risk is essential for timely intervention and improved patient outcomes. Our study aims to develop a predictive model for quantifying risk of irAEs before their clinical diagnosis, utilizing biomarkers derived from medical imaging. Such an approach should provide a framework for continuous risk assessment throughout cancer treatment. To achieve this, we employ logistic regression and Bayes theorem to model risk probabilities, focusing on the transition between biomarker values that predict irAE onset and those that do not.
%
Using quantitative imaging features extracted from 18F-FDG PET/CT images of 58 metastatic melanoma (MM) patients treated with anti-PD-1 and/or anti-CTLA-4 ICIs, we identified biomarker values as independent variables in our risk prediction models. Our analysis targets three organs frequently affected by irAEs: the bowel, lungs, and thyroid. For each method, we calculated optimal and median models, along with $95\%$ confidence intervals, using maximum likelihood estimation (MLE) and bootstrapping to address the small sample size.
%
Our results demonstrate that both models produce meaningful risk quantification, visualized as S-shaped curves, with an inflection point representing the critical threshold between low and high irAE risk. The optimal models (for each organ separately) meaningfully describe the data within the range of biomarker values. The methods presented offer a robust means of quantifying irAE risk in real-time, enabling proactive patient monitoring and by integrating such predictive tools into clinical practice, we aim to enhance early detection of irAEs and improve treatment outcomes for cancer patients undergoing ICI therapy. This approach holds the potential to improve clinical outcomes significantly.
\end{abstract}

\noindent{\it Keywords\/:} risk quantification, immune-related adverse events (irAE), biomarker value, SUV percentile, modelling, toxicity, uncertainty

\submitto{\PMB}
\maketitle

\section{Introduction}
This research is focused on developing a methodology that quantifies the risk of the therapeutic process with an analysis of the uncertainty of the input data and presented methods. 

Quantitative risk assessment (QRA) is a systematic approach to evaluating risks by assigning numerical values to different potential risks, typically involving probabilities of events and their associated impacts. It’s used in fields like finance, engineering, environmental science, and health and safety to assess the likelihood and consequences of risk scenarios and make informed decisions. QRA can include different techniques, such as mathematical modeling, numerical data, and statistical analysis, that can evaluate the probability and impact of different factors on the outcome. The development of risk models dates back to the mid-20th century when their application was more prominent in other disciplines. Several measures for evaluating risk model performance have been developed across fields such as meteorology, economics, and psychology \cite{brier1950verification, murphy1973new, yates1982external, wilks2011statistical, toma2014quantitative}. Our focus is primarily on the application and implementation of developed models in medicine \cite{pepe2003statistical, zhou2014statistical, steyerberg, bellazzi2008predictive}. However, this was not so common in branches of medicine, specifically not in cancer treatment analysis, aimed at the risks of treatment.  The most developed example of the use of quantitative information for risk assessment in medicine is drug approval practice, where different frameworks have been developed to improve transparency in decision-making \cite{guo2010review, kurzinger2020structured, lineberry2016recommendations}. Also, several studies have been conducted on risk prediction models aimed at determining whether a patient is likely to develop cancer or experience a recurrence in the future \cite{richter2018review, kim2012development, yu2016development}. In such a case, the risk model is a statistical approach used to estimate an individual's probability of experiencing a future adverse outcome within a specified time frame, based on clinical biomarkers, demographic, and lifestyle information. These models are becoming increasingly significant in medical practice as clinical care shifts towards greater personalization based on individual characteristics and needs. Recent literature provides valuable insights into the development, application, and evaluation of such risk models \cite{janes2008assessing, cook2008statistical, pepe2010potential, huang2009parametric}. Modeling disease risk or recurrence is naturally approached as a survival analysis problem, with many studies employing survival analysis techniques to develop predictive models. The Cox Proportional Hazards model \cite{cox2018analysis} is commonly used for this purpose, as it accommodates time censoring and supports multivariate analysis. 
%This regression model generates a time-dependent function based on baseline covariate values, estimating the probability of an event occurring at any future point in time.
Logistic regression (LR) is another widely used statistical model, particularly suited for binary outcome prediction \cite{khoshgoftaar1999logistic, nick2007logistic}. This method supports multivariate analysis by building a linear regression model based on the covariates and applying a logistic function to distinguish between the two output classes. It is particularly appropriate for models involving disease states (diseased or healthy) and decision-making (yes or no), and therefore is widely used in studies in the health sciences \cite{shipe2019developing, boateng2019review, schober2021logistic}. Comparing the occurrence of adverse events based on a covariate of interest (e.g., treatment) is a common question in drug safety analysis \cite{pillans2008clinical, murff2003detecting}. Regarding the statistical aspects of adverse event (AE) analysis, guidelines rarely address this topic in depth, instead primarily concentrating on the collection and reporting of data \cite{lineberry2016recommendations}. A recent publication by Coz et al. \cite{coz2024overview} provides an overview of existing regression models suitable for comparing adverse events and discusses the selection of models in relation to the specific characteristics of the adverse events under consideration, where logistic regression and Hazards model stand out as the most prevalent.

Uncertainty estimation is essential for generating confidence evaluations alongside model predictions. The nature of uncertainties and the approaches to addressing them have long been subjects of discussion among statisticians, scientists, engineers, and other specialists \cite{benjamin2014probability, faber2005treatment, lindley2000philosophy, pate1996uncertainties}. This is especially significant in medical imaging, where assessing uncertainty in the model's predictions can help identify areas of concern or offer supplementary information to clinicians \cite{zou2023review}. Uncertainty estimation provides a confidence score, enabling users to assess the reliability of the model's output and recognize instances where the model's performance may be suboptimal \cite{abdar2021review}. Incorporating uncertainty information improves the decision-making process, mitigates risk, and ensures the appropriate involvement of medical professionals in critical cases. There are two main types of uncertainty to be quantified: data uncertainty (aleatoric), which comes from noise in data, including measurement inaccuracies or labeling errors, and model uncertainty (epistemic), which comes from a lack of knowledge or information about the underlying model or data distribution, or from an inadequate model structure \cite{zou2023review, der2009aleatory}. Most of the problems include both types of uncertainties. Aleatoric uncertainty can be quantified by training the model to produce a distribution of possible predictions instead of a single-point estimate, while epistemic uncertainty can (theoretically) be minimized by employing more complex models, gathering additional data, or applying regularization techniques \cite{zou2023review, der2009aleatory, field2007model, kiureghian1989measures}.

Immune Checkpoint Inhibitors (ICI) have significantly improved outcomes and the median overall survival (OS) rate in patients with metastatic melanoma and a variety of other malignancies \cite{wolchok2017overall}, with the cost of high risk of severe serious side effects, named immune-related adverse events (irAE), \cite{gandy2020immunotherapy, lang2019clinical, eshghi201818f, wang2018fatal}. Early detection of irAE is critical to enable optimal treatment planning, minimize treatment interventions, avoid unwanted complications, and improve clinical outcomes. Disease status and treatment responses to ICI are monitored with whole-body 18F-fluorodeoxyglucose positron emission tomography/computed tomography (18F-FDG PET/CT) \cite{aide2022PETmelanoma,filippiRolePET}, which is also sensitive to inflammation of the pathogenesis of irAE \cite{aide2022PETmelanoma}. 
%Several studies have tried to identify imaging or clinical biomarkers extracted from lesions to assess patient's response to immunotherapy \cite{anwar2018, basler2020, flaus2021, dirks2023, peisen2024os, hindie2020, kudura2021} and correlation of organ uptake and response \cite{sachpekidis2023, iravani2020}.
For the assessment of irAE, clinical biomarkers for the prediction of irAE have been explored in \cite{nadaraja2024}. Besides that, the study by Hribernik et al. \cite{Hribernik2021} appears to investigate imaging biomarkers of immune-related adverse events (irAE). They have developed percentiles of organ FDG uptake distribution as possible biomarkers of irAE, to identify, monitor and manage irAEs in patients undergoing immunotherapy, which have been found to be robust by Huff et al.\cite{Huff2021}. 
%To the best of our knowledge, no other study has looked at imaging biomarkers as predictors of irAE. 

%Currently, there are no methodologies to make timely predictions of the occurrence of clinically significant adverse events and assessment of patient-specific treatment risks.  
The process of immunotherapy treatment requires constant supervision and timely decision-making to ensure that the therapy risk for the patient will not significantly affect the quality of his life. Most studies consider risk by assessing different adverse events arising from therapy \cite{wang2018fatal, puzanov2017managing}. However, there is a lack of quantification of this risk, encompassing mathematical interpretation, and the ability to stop or modify therapy in a timely manner when potentially life-threatening circumstances arise. 
%A model for this purpose could provide valuable insights by establishing thresholds or risk levels that indicate when treatment modifications or discontinuations are necessary to prevent severe impacts on patient quality of life. 
%Furthermore, such an approach would support decision-making by enabling continuous reassessment based on patient response, thereby allowing the therapy to adapt dynamically to individual needs. 
%Consequently, there is a critical need for risk quantification to facilitate timely responses to therapy adverse events.

In this paper, we aim to propose modeling the probability of developing irAE (risk/toxicity) based on the values of biomarkers that were extracted from patient's images before and during the therapeutic process. The main idea behind this is to develop a model that could predict the appearance of side effects of therapy (irAE) before they are clinically diagnosed, which will support doctors and patients during the therapeutic treatment. Using data from a retrospective study, we make a predictor model for probabilities based on imaging findings (biomarkers) and evaluate its uncertainties. This model should make a reasonable bound between biomarker values that correspond to healthy tissue and the ones that predict the irAE. While our methodologies will be generally applicable in risk analysis, we will focus on the assessment of immunotherapy risks in metastatic melanoma, as a case study. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Method}



\subsection{Patient data}

We analyzed retrospectively collected patient data and $^{18}$F-FDG PET/CT imaging of patients with metastatic melanoma (MM) who were treated per standard of care with ICI (anti-CTLA-4 or/and anti-PD1) at two institutions: the Institute of Oncology Ljubljana (OIL), Slovenia or at the University of Wisconsin Carbone Cancer Centre (UW), Madison, WI, USA. The date of clinical diagnosis and grade of irAE were acquired via chart review.  

Images of 58 MM patients were retrospectively analyzed for modeling the risk assessment, by following the occurrence of irAE and biomarker selection. We analyzed three target organs, which are most commonly affected by irAE: the lung, bowel, and thyroid. All patients were divided into two groups for each organ of interest; they were assigned to the NC (Normal Control) group, if they did not experience irAE in the organ of interest, or to the AE (Adverse Event) group if they experienced an adverse event in the organ of interest. More about patient data and information on image acquisition can be found in the past work by Hribernik et al. \cite{Hribernik2021}, where the same patient population was used.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{[18]F-FDG PET/CT image acquisition and analysis}

In this patient cohort, 18F-FDG PET/CT were primarily acquired for treatment response assessment. Images were acquired on five PET/CT scanners; after reconstruction, images were normalized by patient weight and injected dose to compute standardized uptake values (SUV). 
%More information on image acquisition can be found in \cite{Hribernik2021}.
To detect adverse event in an organ of interest, Quantitative Imagining Biomarkers (QIB) extracted from patient's $^{18}$F-FDG PET/CT need to be developed. In our study, we utilized percentiles of organ FDG uptake distribution (SUV percentiles) as potential biomarkers for identifying, monitoring, and managing irAEs in patients undergoing immunotherapy, as these have been demonstrated to be effective descriptors in previously published studies \cite{Hribernik2021, Huff2021}. 

The organs of interest in this paper (bowel, thyroid, lung) were segmented using a Convolutional Neural network (CNN) to enable $^{18}$F-FDG organ uptake quantification \cite{Hribernik2021}. Percentiles of standardized uptake value  ($SUV_\%$) distribution were extracted from each $^{18}$F-FDG PET for all three organs of interest. We focused on the following percentiles, that were found optimal to detect irAE in Hribernik et al. \cite{Hribernik2021}: $SUV_{95\%}$ for bowel and lungs, and $SUV_{75\%}$ for thyroid. For patients who had multiple $^{18}$F-FDG PET/CT scans, the maximum value of ($SUV_\%$) for each organ of interest was used as a predictor of irAE. To model the probability prediction of a patient experiencing an adverse event, $SUV_\%$ in three organs of interest, and their associated clinical data, were analyzed independently. The number of patients in each group sorted by organs of interest is shown in Table \ref{tab:data_sum}.

\begin{table}[!htb]

\centering
\begin{tabular}{c|c|c}
Organ       &   $\nc$ group ($N_\nc$) &  $\ae$ group ($N_\ae$)   \\ 
\hline
bowel       &  52 & 6  \\
lung        &  53 & 5  \\
thyroid     &  49 & 9
\end{tabular}

\caption{Patient population ($N=58$) used in the study, separated into normal control ($\nc$) and adverse event ($\ae$) group for each analysed organ,  $N=N_\nc + N_\ae$.}
\label{tab:data_sum}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Mathematical modelling probability of toxicity and its uncertainty }\label{sec:matmodel}

We consider a random variable $(X,Y)$, where $X$ represents a biomarker value and $Y \in \{\nc, \ae\}$ its binary associated organ's state: normal control and adverse event case.  In the process of measurements, we obtained a sample of this random variable:
%
\begin{equation}
  {\cal S} = \left\{(x_i,y_i) : y_i \in \{\ae,\nc\},~i = 1,\ldots,N\right\}\>,
\end{equation}
%
where $N$ represents number of patients (in our case $N=58$), $x_i$ is the biomarker value ($SUV_\%$) for the $i$th patient and the considered organ, while $y_i$ is the organ state (NC or AE).
In our model, prior information on which we based conditional probability is the value of biomarker ($SUV_{\%}$), denoted by $X$. We are interested in the probability of adverse events based on the value of the selected biomarker obtained from $^{18}$F-FDG PET/CT images. These biomarker values should be specific to the observed disease and differ in the case of diseased and healthy tissue, \cite{Hribernik2021, Huff2021}. As our study includes irAE in three organs, the lung, bowel, and thyroid, for each of them a specific $SUV_{\%}$ was used, found optimal to detect irAE in Hribernik et al. \cite{Hribernik2021}.

The conditional probability of encountering adverse events as a function of biomarker value, can be denoted as $P(Y = \ae \mid X = x)$.

Let us denote a particular model for conditional probability as $f(x \mid \boldsymbol{\beta})$, where $\boldsymbol{\beta}$ is a vector of parameters. The optimal model or its associated parameters $\boldsymbol{\hat\beta}$ are found using the maximum likelihood estimation (MLE) approach. Due to sample finiteness, the model parameters $\boldsymbol{\beta}$ can be determined up to their probability distribution, because we are considering only a sample from the population \cite{Efron1982}. The latter is obtained by a combination of parametric and non-parametric bootstrapping techniques \cite{Efron1982}. 

We see that values of model parameters $\boldsymbol{\beta}$ is a random vector and consequently, the model value $f(x \mid \boldsymbol{\beta})$ at given $x$ is a random variable. The uncertainty of model values is represented by their confidence intervals (CIs). The CI of model values $[L_\alpha(x), U_\alpha(x)]$ with significance level $\alpha$ at a specific biomarker value $x$ is determined by the condition:
%
\begin{equation}
    {\rm Prob}\left[L_\alpha(x) \le f(x \mid \boldsymbol{\beta}) \le U_\alpha(x)\right] 
    = 1 - \alpha \>,
     \label{eq:conf_int}
\end{equation}
%
where $ {\rm Prob}[]$ denotes the probability of events satisfying the condition stated in brackets. Here we focus only on equal-tail CI that fulfil additional constraints:
%
\begin{equation}
    {\rm Prob}\left[f(x \mid \boldsymbol{\beta}) \le L_\alpha(x)\right] = 
    {\rm Prob}\left[f(x \mid \boldsymbol{\beta}) \ge U_\alpha(x)\right] = \frac{\alpha}{2} \>.
\end{equation}
%
Another interesting quantity is the median of model values, denoted by $M(x)$, and is determined by the condition:
%
\begin{equation}
    {\rm Prob}[f(x \mid \boldsymbol{\beta}) \le M(x) ] = \frac{1}{2} \>.
\end{equation}
%
Median $M(x)$ is typically near to optimal model $f(x \mid \boldsymbol{\hat\beta})$, but diverges, where $f(x \mid \boldsymbol{\hat\beta})$ is near $0$ or $1$, i.e. boundaries of the model value domain. CIs are numerically determined using parametric (estimating distribution of parameters) and non-parametric (sample with replacements) bootstrapping techniques via basics percentile method, see e.g. \cite[Ch.~5.3.1]{DavisonHinkley1997}.

The bootstrapping technique often estimates a population's properties by resampling with replacing the original data and fitting the model to it. So each new sample is associated with corresponding model parameters, and we end up with a set of model parameters. In the non-parametric bootstrapping approach, the set of model parameters is used directly to calculate model values at conserving values of independent variables, and the confidence interval (CI) of model values is obtained as quantiles of the model values, i.e. basics percentile method. On the other hand, in the parametric bootstrapping approach, we approximate an analytical distribution of parameters for the set of model parameters and use a random generator for that analytical distribution to generate parameters and corresponding model values.  Based on this, we calculate the CI using the basic percentile method.

We considered two fundamentally different mathematical models, i.e. logistic/logit regression and binary group Bayesian modeling. In the latter, within Bayesian modeling, the groups (NC and AE) are modeled using normal and log-normal distributions. For these model types, we report on optimal model $f(x \mid \boldsymbol{\hat\beta})$, model median $M(x)$, and confidence intervals/bands ${\rm CI}(x)$.

\subsubsection{Logistic regression modeling} \label{sec:log_met}

Logistic regression is a statistical method used to model the relationship between independent variables and a binary outcome. In this study, the binary outcome is defined as the presence or absence of adverse events. The logistic regression framework establishes a linear relationship between the independent variables and the natural logarithm of the odds of the outcome, which corresponds to a sigmoid relationship between the independent variables and the probability of the outcome being 1. This ensures that the predicted probabilities remain constrained within the interval [0, 1].
The model assumes that the odds of experiencing adverse events are a linear function of the biomarker values, defined by the intercept and rate parameters.

Our objective is to determine the conditional probability $P(Y = \ae  \mid X = x)$ using logistic regression performed via the MLE method. To address this, we make use of a logistic function, defined as:
%
\begin{equation}
  f(x  \mid \boldsymbol{\beta}) = \frac{1}{1 + e^{-\beta_0 - \beta_1 x}}\>,
  \label{eq:logistic_fun}
\end{equation}
%
with parameter vector $\boldsymbol{\beta} = [\beta_0, \beta_1]^\mathsf{T}$, where $\beta_0$ and $\beta_1$ represent the intercept and rate parameter, respectively, and the conditional probability of irAE occurrence is $P(Y = \ae  \mid X = x)=f(x  \mid \boldsymbol{\beta})$.

The optimal parameter $\boldsymbol{\beta}$ is obtained using MLE method, i.e., by  maximizing the log-likelihood function of the given problem:
%
\begin{equation}
    \log L(\boldsymbol{\beta}) 
    = \sum_i
        \left[  \delta_{y_i, \ae}\log(f(x_i \mid \boldsymbol{\beta})) + 
                (1-\delta_{y_i, \ae})\log(1-f(x_i \mid \boldsymbol{\beta}))
        \right] \>,
    \label{eq:logistic_llnf}
\end{equation}
%
where $\delta_{i,j} = (1:i = j;~0:\textrm{otherwise})$ is Kronecker delta. The local optimum $\boldsymbol{\hat \beta}$ satisfies the zero gradient equation 
%
\begin{equation}
    \nabla_{ \boldsymbol{\hat\beta}} 
    \log L(\boldsymbol{\hat\beta}) 
    = \sum_i \left[\delta_{y_i, \ae} - f(x_i \mid \boldsymbol{\hat\beta})\right] 
    \begin{bmatrix} 1 \\ x_i \end{bmatrix}
    = \mathbf{0} \>,
\end{equation}
%
and has a Hessian matrix given by
%
\begin{equation}
    \begin{split}
    {\bf H}(\boldsymbol{\hat\beta}) 
        &= \nabla_{\boldsymbol{\hat \beta}} \nabla_{\boldsymbol{\hat \beta}}^\mathsf{T}  \log L( \boldsymbol{\hat\beta}) \>,\\
        &= -\sum_i \left[1 - f(x_i \mid \boldsymbol{\hat\beta})\right] 
                   f(x_i \mid  \boldsymbol{\hat \beta}) 
        \begin{bmatrix} 
            1 & x_i  \\
            x_i & x_i^2  
        \end{bmatrix} \>.
    \end{split}
    \label{eq:logistic_hessian}
\end{equation}
%
In our study, we find the global optimal parameters $\boldsymbol{\hat \beta}$ using L-BFGS algorithm \cite{Nocedal2006} with a meaningful initial guess. The inverse of the Hessian matrix ${\bf H}(\boldsymbol{\hat\beta})$, Eq.~\eqref{eq:logistic_hessian}, is an asymptotic estimate of the variance-covariance matrix of parameters $\boldsymbol{\Sigma}_{\rm log}$ in MLE method (for details, see \cite{Newey1994} and \cite[Chapters 7.2 and 7.3]{Lehmann2004}). This means that for large sample sizes, we expect that the parameters are approximately normally distributed as
%
\begin{equation}
    \boldsymbol{\beta} \sim N(\boldsymbol{\hat\beta}, \boldsymbol{\hat\Sigma}_{\rm log}),\qquad
    \boldsymbol{\hat\Sigma}_{\rm log} = -{\bf H}(\boldsymbol{\hat\beta})^{-1} \>,
    \label{eq:mvn_log}
\end{equation}
%
where $N(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ represents a multivariate normal distribution with mean $\boldsymbol{\mu}$ and variance-covariance matrix $\boldsymbol{\Sigma}$. %
%

The uncertainty estimation is determined through model confidence intervals, which are constructed using a bootstrapping technique, with several known variants available \cite{Efron1982}. This was done by applying the bootstrapping technique to estimate the parameters of the distribution, which subsequently resulted in uncertainties within the resulting model estimation (see \cite{Efron1994} and \cite[Chapter 15.6]{Press1992}).  Due to the small sample size of group $\ae$, non-parametric bootstrapping (resampling with replacement) produces a parameter distribution with multiple competing peaks. Consequently, parametric bootstrapping is required, as described in detail below.

In the case of logistic regression modeling, the non-parametric bootstrapping results in a distribution of parameters with several peaks, which are similarly probable and one dominant. This can be seen in Figure \ref{fig:logit_boots}.a and its zoom into dominant peak \ref{fig:logit_boots}.b. Each peak corresponds to a different number of AE patients in the samples obtained via resampling with replacement technique. The dominant peak is associated with the case described by original data, and at approximately its position we find the MLE parameter $\boldsymbol{\hat\beta}$ marked by a big red dot. The red lines in side figured of \ref{fig:logit_boots}.b represent the marginalized probability density of the parameters Eq.~\eqref{eq:mvn_log}, which only captures the general position and roughly the shape of the bootstrapped peak.
%
\begin{figure}[!htb]
\centering
\includegraphics[width=7cm]{logit_boots}
\includegraphics[width=7cm]{logit_boots_sel}
\caption{Distribution of logistic regression parameters for lung data obtain via non-parametric bootstrapping (a) and zoom into the region of the highest peak (b). In the latter, the red dot marks the MLE parameter, and the red lines in the side figures represent the marginalized asymptotic normal multivariate distribution of the parameters.}
\label{fig:logit_boots}
\end{figure}
%
Because the distribution of parameters obtained by bootstrapping has a lot of finite sample size anomalies, which are difficult to interpret and should vanish in large sample size, we instead use asymptotic limit results and approximate the distribution of parameters $f_{\boldsymbol{B}}(\boldsymbol{\beta})$ as a multivariate normal distribution given by Eq.~\eqref{eq:mvn_log}. This is then used in parametric bootstrapping: sampling parameters from approximated distribution and observe resulting model values. 


\subsubsection{Bayesian probability modeling}

The other approach includes Bayesian modeling to access the required probability, which makes this model more complicated than the logistic regression modeling. Bayesian modeling, with its formal use of prior information and iterative updates based on accumulating knowledge, naturally supports decision theory.  Bayes' theorem enables modeling depending on previous outcomes and thus provides complete information in such types of modeling of patients’ conditions, the so-called conditional probabilities whose realizations are dependent and induced by previous events. Bayesian modeling uses Bayes' theorem to describe the conditional probability of an event based on data as well as prior information or beliefs about the event or different conditions that are related to the event. 

The joint probability distribution of the pair of variables $(X, Y)$ is
%
\begin{equation}
  p_{X,Y}(x,y) =  p_X(x \mid y) P_Y(y) \>,
\end{equation}
%
where $p_X(x \mid y)$ is the conditional probability density function, and $P_Y(y)$ is the discrete value, respectively represents the probability that the patient is in one of the observed groups $y  \in \{\ae, \nc\} $, based on data from a retrospective study.
The total probability law gives the probability density function of $X$:
%
\begin{equation}
 p_X (x) = p_X(x \mid \ae) P_Y(\ae) + p_X(x \mid \nc) P_Y(\nc) \>.
\label{eq:tot_prob}
\end{equation}

Namely, we are interested in the conditional probability that a person will experience irAE, i.e. $Y=\ae$, if we know the biomarker $X$ value observed from PET images ($SUV_\%$ specific for each organ of interest). By using the Bayesian formula, this could be modeled as: 
%
\begin{equation}
  P(Y= \ae  \mid X = x) =  
  \frac{P_Y(\ae) p_X(x \mid \ae)}{ p_X(x \mid \ae) P_Y(\ae) + p_X(x \mid \nc) P_Y(\nc)}\> ,
  \label{eq:cond_prob1}
\end{equation}
%
where the denominator equals total probability, Eq.~\eqref{eq:tot_prob}. By introducing ratio between probability densities $r(x)$ and the odds for a patient to be in $\nc$ group $O$ given by
%
\begin{equation}
    r(x) = \frac{p_X(x \mid \nc)}{p_X(x \mid \ae)} \quad \textrm{and}\quad
    O = \frac{P_Y(\nc)}{P_Y(\ae)} \>,
\end{equation}
%
the conditional probability can be rewritten compactly to
%
\begin{equation}
  P(Y= \ae  \mid X = x) = \frac{1}{1 + O\, r(x) } \>.
  \label{eq:cond_prob2}
\end{equation}
%
This form is particularly convenient to simplify further analysis. 

The challenge here is to model the probability density function of SUV percentiles for both groups, $p_X(x \mid \ae)$ and $p_X(x \mid \nc)$, with the proper distribution.  
The probabilities for a patient being in each of the groups, i.e. $P_Y(\nc)$, $P_Y(\ae)$, are fixed for our study case and are calculated 
using the data in Table~\ref{tab:data_sum} as 
%
\begin{equation}
    P_Y(\nc) = \frac{N_\nc}{N} \quad\textrm{and}\quad P_Y(\ae) = \frac{N_\ae}{N} \>. 
\end{equation}
%
%$N_\nc$ and $N_\ae$ are the numbers of patients in the normal control and adverse effect group, respectively, and $N=N_\nc+N_\ae$ is the number of all patients. 

%[{\bf Katja: How do the values in the table compare to other cases reported in the literature}]


\paragraph{Gaussian type of distributions.} Let us model the conditional probability densities associated with both groups, i.e., $p_X(x \mid \ae)$ and $p_X(x \mid \nc)$, as Gaussian distributions:
%
\begin{equation}
  p_X(x \mid y) = 
    \frac{1}{\sqrt{2 \pi} \sigma_y} 
    \exp\left(-\frac{(x - \mu_y)^2}{2 \sigma_y^2}\right) \qquad
    y  \in \{\ae, \nc\}   \>,
\end{equation}
%
where $\mu_y$ and $\sigma_y$ are the mean and standard deviation of the observed group $y  \in \{\ae, \nc\} $. Then we can write the requested conditional probability as
%
\begin{equation}
  P(Y = \ae  \mid X = x) = \frac{1}{1  + \exp\left(F_{\rm norm}(x)\right)} \>,
  \label{eq:prob_gauss}
\end{equation}
% 
where $F_{\rm norm}$ is given by
%
\begin{align}
  F_{\rm norm}(x) &= \ln\frac{P_Y(\nc)\sigma_\ae}{P_Y(\ae)\sigma_\nc} 
                    -\frac{(x - \mu_\nc)^2}{2\sigma_\nc^2} + 
                     \frac{(x - \mu_\ae)^2}{2\sigma_\ae^2} \>,
   \label{eq:Fgauss}
\end{align}
%
Function $F_{\rm norm}$ is a quadratic polynomial in variable $x$ and so by writing $F_{\rm norm}(x) = a_0 + a_1 x + a_2 x^2$ its coefficients are 
%
\begin{align}
  a_0 &= \ln \left [\frac{P_Y(\nc)\sigma_\ae}{P_Y(\ae)\sigma_\nc} \right]+ 
        \frac{1}{2} \left [-\frac{\mu_\nc^2}{\sigma_\nc^2} + \frac{\mu_\ae^2}{\sigma_\ae^2}\right] \>, \\
  a_1 &=  \frac{\mu_\nc}{\sigma_\nc^2} - \frac{\mu_\ae}{\sigma_\ae^2}\>, \\
  a_2 &=  \frac{1}{2}\left [-\frac{1}{\sigma_\nc^2} + \frac{1}{\sigma_\ae^2}\right] \>.
  \label{eq:gauss_coeff}
\end{align}
%
Notice that the conditional probability has only three independent parameters, but the total number of parameters is four. The formula for conditional probability, Eq.~\eqref{eq:prob_gauss}, has a similar form as the logistic function, Eq.~\eqref{eq:logistic_fun}, but has significantly different properties. Generally, standard deviations $\sigma_\ae$ and $\sigma_\nc$ are different and in this case asymptotic behavior is determined by the sign of coefficients $a_2$: 
%
\begin{equation}
 \lim_{x\to\pm \infty} P(Y = \ae  \mid X = x) = 
 \left \{
     \begin{array}{lll}
     1 & :& \sigma_\nc < \sigma_\ae \\
     0 &: & \sigma_\nc > \sigma_\ae
     \end{array} \right. \>.
  \label{eq:bayes_normal_asymp}
\end{equation}
%
This means generally conditional probability $P(Y=\ae  \mid X = x)$ has either 1 or 0 asymptote and so it is not an S-shaped curve. If the standard deviations are identical, i.e.  $\sigma_\nc = \sigma_\ae \equiv \sigma$, $F(x)$ becomes linear function and conditional probability $P(Y=\ae  \mid X = x)$ is translated to a logistic model:
%
\begin{equation}
   P(Y = \ae  \mid X = x) = \frac{1}{1 + O \exp[b (x - a)]}\>,
\end{equation}
%
with parameters $a = (\mu_\nc + \mu_\ae)/2$, $b = (\mu_\nc - \mu_\ae)/\sigma^2$. In order to have an expected behavior of a monotonically increasing S-shaped curve from 0 to 1 values with increasing value of $x$, parameter $b$ should be negative. %


\paragraph{Log-normal type of distributions.}

In probability theory, a log-normal distribution is a continuous probability distribution of a random variable whose logarithm is normally distributed. If the random variable $X$, in our case SUV percentile, is log-normally distributed, then $ln(X)$ has a normal distribution, \cite{scarpelli2016we,thie2000diagnostic}. Mainly, normal distributions can allow for negative random variables while log-normal distributions include all positive variables. The values of the mentioned biomarker that we are considering through this paper (SUV percentile) can only be positive, which implies that this fact can justify the use of the log-normal distribution for modeling the probability density function of the conditional SUV percentile distribution ($p_X(x \mid \ae), p_X(x \mid \nc)$). The probability density function for log-normal is slightly different from the Gaussian.

 If we model the conditional probability densities associated with both groups, i.e., $p_X(x \mid \ae)$ and $p_X(x \mid \nc)$, as log-normal distributions:
%
\begin{equation}
  p_X(x \mid y) = 
    \frac{1}{x\sqrt{2 \pi} \sigma_y} 
    \exp\left(-\frac{(\ln x - \mu_y)^2}{2 \sigma_y^2}\right) \qquad
    y  \in \{\ae, \nc\}   \>,
\end{equation}
%
where $\mu_y$ and $\sigma_y$ are the mean and standard deviation of the observed group $y  \in \{\ae, \nc\} $, then we can write the requested conditional probability as
%
\begin{equation}
  P(Y = \ae  \mid X = x) = \frac{1}{1  + \exp\left(F_{\rm norm}(\ln x)\right)} \>,
  \label{eq:prob_log_normal}
\end{equation}
% 
where $F_{\rm norm}(x),$ defined in Eq. \eqref{eq:Fgauss}, represents a quadratic polynomial derived from the normal distribution to describe groups. In this case, the conditional probabilities are identical to those obtained in the normal distribution scenario, except that the argument is replaced by its natural logarithm, $\ln x$. This correspondence implies that the conditional probability exhibits similar asymptotic behavior:
%
\begin{equation}
 \lim_{x\to 0, \infty} P(Y = \ae  \mid X = x) = 
 \left \{
     \begin{array}{lll}
     1 & :& \sigma_\nc < \sigma_\ae \\
     0 &: & \sigma_\nc > \sigma_\ae
     \end{array} \right. \>.
  \label{eq:bayes_lognor_asymp}
\end{equation}
%
Additionally, similar to the previous case, the conditional probability is generally not an S-shaped curve. It takes an S-shape only when the spreads are equal, i.e. $\sigma_\nc = \sigma_\ae \equiv \sigma$:
%
$$
    P(Y = \ae  \mid X = x)  = \frac{1}{1 + O \exp(-a b) x^b}\>.
$$
%
The derived formula is functionally equivalent to a logistic model in which the logarithm $\ln x$ serves as the argument. The conditional probability exhibits an S-shaped curve that transitions from 0 to 1 as the argument increases, only for negative values of $b$.

\subsubsection{Properties of the Bayesian approach} 
\label{sec:Properties of the Bayesian approach}
In our Bayesian probability modeling, we adopt a two-step process to estimate the conditional probability distribution:
%
\begin{enumerate}
    \item Group-wise fitting: We first fit a probability distribution to each group separately, obtaining parameter estimates for each group.
    \item Combining parameters: The obtained parameters are then combined to express the conditional probability distribution.
\end{enumerate}
%
To approximate the distribution of the parameters of the conditional probability function, we extend this approach by employing a bootstrapping strategy:

\begin{enumerate}
    \item Non-parametric bootstrapping: We perform non-parametric bootstrapping to repeatedly sample data within each group and refit the distributions to obtain sets of bootstrapped parameters for each group.
    \item Parameter combination: We then take the Cartesian product (or equivalently, the outer product with a join operation) of the parameter sets from each group. This step results in a combined set of parameters that characterize the conditional distribution.
\end{enumerate}
%
This approach strongly resembles stratified non-parametric bootstrapping, with the difference that stratification is applied during the fitting of the distributions for each group and not during the fitting of the conditional distribution.

It is important to note that any Bayesian model based on normal or log-normal distributions for describing both groups results in conditional probabilities defined by three independent parameters. Without additional constraints on these parameters, the resulting probabilities are not S-shaped on the domain argument. Imposing these constraints transforms the conditional probability model into a logistic model.

To achieve an $S$-shaped conditional probability that transitions from 0 to 1 as the argument increases, it is often necessary to use less commonly employed distributions, which typically have a greater number of parameters. These distributions can be challenging to justify, especially in scenarios with limited data, as the additional parameters may not be supported statistically. Below, we provide two simple examples of distributions defined on the domain $[0, \infty)$:

\begin{enumerate}
    \item Beta-prime (inverted Beta) distribution \cite{Johnson1995}: Consider the beta-prime distribution with the probability density function:
    %
    \begin{equation}
        f_{\text{beta'}}(x \mid \alpha, \beta) = 
        \frac{x^{\alpha - 1} (1 + x)^{-\alpha - \beta}}{\mathrm{B}(\alpha, \beta)} 
        \qquad \alpha, \beta > 0,
    \end{equation}
    %
    where $\mathrm{B}(\alpha, \beta)$ is the beta function. This distribution can be used to express the probability densities for the two groups as follows:
    %
    \begin{equation}
        p_X(x \mid y) = f_{\text{beta'}}(x \mid \alpha_y, \beta_y) \qquad
        y \in \{\ae, \nc\} \>.
    \end{equation}
    %
    The conditional probability $P(Y = \ae \mid X = x)$ takes an $S$-shaped form if the parameters satisfy the following conditions:
    %
    \begin{equation}
        \alpha_\nc < \alpha_\ae \quad \text{and} \quad \beta_\nc > \beta_\ae.
    \end{equation}

    \item Gamma distribution \cite{Johnson1994}: Next, consider the gamma distribution with the probability density function:
    %
    \begin{equation}
        f_{\text{gamma}}(x \mid \alpha, \lambda) = 
        \frac{\lambda^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\lambda x} 
        \qquad \alpha, \lambda > 0,
    \end{equation}
    %
    where $\Gamma(\alpha)$ is the gamma function. This distribution can be used to describe the probability densities for the two groups as:
    %
    \begin{equation}
        p_X(x \mid y) = f_{\text{gamma}}(x \mid \alpha_y, \lambda_y) \qquad
        y \in \{\ae, \nc\} \>.
    \end{equation}
    %
    The conditional probability $P(Y = \ae \mid X = x)$ is $S$-shaped if the parameters satisfy the conditions:
    %
    \begin{equation}
        \alpha_\nc < \alpha_\ae \quad \text{and} \quad \lambda_\nc > \lambda_\ae.
    \end{equation}
\end{enumerate}
%
In the two mentioned distributions, additional parameters—location and scale—are also involved. As a result, the modeling would depend on four parameters, which, given our problem and the available data size, is impractical.
For further information about univariate distributions in use, see e.g. \cite{Johnson1995}  and \cite{Johnson1994}.

\section{Results}

%We present two different approaches, logistic and Bayesian, to model the probability of encountering adverse effects as a function of biomarker value, $SUV_\%$.  
We are focused on two problems: obtaining the optimal model in the MLE sense and estimating the confidence intervals of the model values. 

\subsection{Logistic regression modeling}
We find the optimal value in the MLE sense, and this gives a sigmoid-shaped relationship between the biomarker value $SUV_\%$ (independent variable) and the probability of risk (dependent variable). The optimal model, median model, and 95\% confidence intervals for our data's model values are computed as explained in Sec.~\ref{sec:log_met}. The results are presented in Figure \ref{fig:logit_mvn} and optimized parameter values for each organ in Table \ref{tab:logit_res}. In all cases, lung, bowel, and thyroid data, the optimal solutions were successfully obtained and matched the median of the distribution of model values. 
%
\begin{table}[!htb]
\sisetup{round-mode=places,round-precision=3}
\caption{Parameters of the optimal logistic regression model.}
\label{tab:logit_res}
\centering
\begin{tabular}{l|l|l}
organ   & $\hat \beta_0$ & $\hat\beta_1$ \\
\hline
lung    & \num{-11.683978} & \num{5.687027} \\
bowel   & \num{-3.474206}  & \num{0.337083} \\
thyroid & \num{-4.387076}  & \num{1.163232}
\end{tabular}
\end{table}

The optimal solution aligns with meaningful interpretations of the logistic model to discriminate between groups based on biomarker values ($SUV_\%$). The width of the confidence intervals (CIs) is influenced by the clarity with which patients can be assigned to groups (NC and AE). Thus, the width of the CIs reflects the degree of separation between the groups within the dataset. Consequently, in the lung dataset, where the $\nc$ and $\ae$ groups are most distinctly separated, the CIs are the narrowest.

The disadvantage of this method is that the classes are not pre-separated; data from both groups (NC and AE) are analyzed collectively, and the model determines the division based on the data while optimizing the parameters. On the one hand, this approach simplifies the model, but on the other hand, it impacts the optimization of the parameters.
%
\begin{figure}[!htb]
\centering
    \includegraphics[width=0.49\textwidth]{logit_mvn_lung}%
    \includegraphics[width=0.49\textwidth]{logit_mvn_bowel}\\
    \includegraphics[width=0.49\textwidth]{logit_mvn_thyroid}
\caption{Results for the conditional probability $P(Y=\ae \mid X=x)$ obtain via logistic regression and estimated model uncertainty -- 95\% two-sided confidence intervals of model values for data connected to lungs (top left), bowel (top right) and thyroid (bottom). }
\label{fig:logit_mvn}
\end{figure}

\subsection{Bayesian probability modeling}

Compared to the previously mentioned logistic modeling, Bayesian modeling represents a more advanced theoretical approach. It incorporates assumptions about the distributions of biomarker values within each group (NC and AE) individually, as well as the parameters defining these distributions. These parameters are subject to constraints, making the optimization process numerically sensitive. The primary challenge lies in modeling these distributions, as the dataset of biomarker values available for analysis is very limited - particularly in the AE group, which contains only 6-9 values, depending on the organ of interest. Consequently, we restricted our analysis to distributions with a maximum of two parameters.



\begin{table}[!htb]
\sisetup{round-mode=places,round-precision=3}
\caption{Parameters of the optimal Bayesian model with normal distribution describing conditional probability of groups.}
\label{tab:bayesian_gauss}
\centering
\begin{tabular}{l|l|l|l|l}
organ   & $\hat \mu_\nc$ & $\hat\sigma_\nc$ & $\hat \mu_\ae$ & $\hat\sigma_\ae$\\
\hline
lung    & \num{1.30761063}  & \num{0.25845619} & 
           \num{2.44364713} & \num{0.88137532} \\
bowel   & \num{3.29019878} & \num{1.61425312} & 
          \num{4.86590897} & \num{2.14788649} \\
thyroid & \num{1.93133207} & \num{0.86266287} & 
          \num{3.07875405} & \num{0.80272582}
\end{tabular}

\end{table}


%
\begin{figure}[!htb]
\centering
    \includegraphics[width=0.49\textwidth]{Bayes_bs_Gauss_lung}%
    \includegraphics[width=0.49\textwidth]{Bayes_bs_Gauss_bowel}\\
    \includegraphics[width=0.49\textwidth]{Bayes_bs_Gauss_thyroid}
\caption{Results for the conditional probability $P(Y=\ae \mid X=x)$ obtain via Bayesian probability modeling and Gaussian type of distribution for both conditional density functions ($p_X(x \mid \ae), p_X(x \mid \nc)$), with estimated model uncertainty by plotting $95\%$ two-sided confidence intervals of model values for data connected to lungs (top left), bowel (top right) and thyroid (bottom). }
\label{fig:Gauss_bs}
\end{figure}

We compute the optimal model, median model, and $95\%$ confidence intervals of the model values for our data connected with three organs of interest. The results are shown in Figure \ref{fig:Gauss_bs} and optimized parameter values in Table \ref{tab:bayesian_gauss}. From the graphs in Figure \ref{fig:Gauss_bs}, it can be seen that for all cases, lung, bowel, and thyroid data, the optimal solutions were obtained, and the model curve is not completely aligned with the median of distribution values. The width of the model value CI is determined by the unambiguity of attributing a patient into groups based on the biomarker value. Therefore, the width of CIs reflects how well the groups are separated in data. Consequently, in the lung dataset, where $\nc$ and $\ae$ groups are the most separated, the CIs are the narrowest in the part of the graph where we have samples, and going further right it becomes wider because there are no values to direct it. The huge reverse back at the beginning of the lung graph is a consequence of quadratic polynomial function $F_{\rm norm}(x)$, Eq.~\eqref{eq:Fgauss}. In the remaining two organs we are considering, the classes overlap a lot, which causes the model line and median to differ significantly. In all three cases, the right side of the graphs is characterized by a wide confidence interval due to the lack of biomarker values that would better determine it.


Figure \ref{fig:Log_normal_bs} represents results for Bayesian modeling and log-normal distribution type of conditional probability density functions. As in previous cases, we compute the optimal model, median model, and 95\% confidence intervals for the model values for our data. Table \ref{tab:bayesian_lognormal} gives the optimized parameter values for each organ of interest. The optimal solution resembles meaningful solutions for discriminating between groups based on biomarker value. In lung data, the optimal solution was successfully obtained and matched the median distribution of model values. In this class, $\nc$ and $\ae$ groups are the most separated, which results in the narrowest CIs. Again, there is a reverse back at the beginning of the graph, for SUV values close to zero, which is attributed to the form of the quadratic function, $F_{\rm norm}(\ln x)$. %Eq. \eqref{eq:Fgauss} with argument $\ln x$. 
The width of CIs reflects how well the groups are separated in data. Consequently, in the bowel and thyroid datasets, where $\nc$ and $\ae$ groups overlap, the CIs are significantly wider. Also, in these groups, especially in bowel data, the median of distribution of model value deviates from the optimal solution. 

\begin{table}[!htb]
\sisetup{round-mode=places,round-precision=3}
\caption{Parameter of the optimal Bayesian model with log-normal distribution estimating conditional probability associated to the groups.}
\label{tab:bayesian_lognormal}
\centering
\begin{tabular}{l|l|l|l|l}
organ   & $\hat \mu_\nc$ & $\hat\sigma_\nc$ & $\hat \mu_\ae$ & $\hat\sigma_\ae$\\
\hline
lung    & \num{0.25012563}  & \num{0.18862822} & 
          \num{0.84495162} & \num{0.3054373} \\
bowel   & \num{1.11554438} & \num{0.35077778} & 
          \num{1.49665276} & \num{0.40295539} \\
thyroid & \num{0.59700658} & \num{0.32538502} & 
          \num{1.0869074} & \num{0.28239909}
\end{tabular}
\end{table}



\begin{figure}[!htb]
\centering
    \includegraphics[width=0.49\textwidth]{Bayes_bs_log_normal_lung}%
    \includegraphics[width=0.49\textwidth]{Bayes_bs_log_normal_bowel}\\
    \includegraphics[width=0.49\textwidth]{Bayes_bs_log_normal_thyroid}
\caption{Results for the conditional probability $P(Y=\ae \mid X=x)$ obtain via Bayesian probability modeling and log-normal type of distribution for both conditional density functions $(p_X(x \mid \ae), p_X(x \mid \nc))$, with estimated model uncertainty by plotting $95\%$ two-sided confidence intervals of model values for data connected to lungs (top left), bowel (top right) and thyroid (bottom). }
\label{fig:Log_normal_bs}
\end{figure}

We naturally aim for $95\%$ confidence intervals (CIs) or bands to be as narrow as possible. However, in our analysis, the CIs across all datasets and models are generally wide. The narrowest CIs occur for biomarker values near the centers of the two groups. Within this range, the median model values align closely with the optimal model. Outside this range, the CIs become significantly wider, rendering them practically ineffective for uncertainty estimation. The overall width of the CIs can be attributed to the inherent characteristics of the data: the AE group is typically very small (6-9 values), resulting in considerable uncertainty in the group discriminators, and there is substantial overlap between the AE and NC groups. The narrowest CIs are observed in the lung dataset, where group separation is the most distinct. In contrast, the CIs are wider in the bowel and thyroid datasets. It is worth noting that, in the logistic regression approach, parametric bootstrapping was employed, with sampled parameters centered around the maximum likelihood estimate (MLE) solution. In contrast, Bayesian modeling utilized non-parametric bootstrapping to derive confidence intervals for the model values.




\section{Discussion}

Our study explores two approaches for modeling the probability of immune-related adverse events based on biomarker values: logistic regression, a common method for binary outcomes, and Bayesian modeling, which can sometimes reduce to logistic regression. We address two key problems: obtaining the optimal model using maximum likelihood estimation (MLE) and quantifying uncertainties through confidence intervals.
There is a notable lack of empirical evaluations for such models and their uncertainty estimation methods, particularly in real clinical scenarios. This limitation makes it challenging to compare and assess different methods and to identify which approaches are most effective and efficient for specific tasks.

\paragraph{Obtaining the optimal model:}
Let's return to our probabilistic model $P(Y = \ae \mid X = x) = f(x \mid \boldsymbol{\beta})$ typically selected by fitting a theoretical distribution to data. While standard goodness-of-fit tests assess overall fit, they often fail to ensure accuracy in the tail, which is critical in structural reliability and risk analysis. Different well-fitted distributions can yield varying probability estimates, introducing epistemic (model) uncertainty. To address this, we parameterize distribution selection, representing model uncertainty through parameter uncertainty. This is achieved using bootstrapping techniques to estimate parameters in our probabilistic approach.

This study examined two theoretically distinct models, yielding results that may seem unexpected at first glance. Although the Bayesian method is more complex and theoretically expected to perform better, logistic regression proved more effective in producing a monotonically increasing S-shaped curve, accurately assigning probabilities to biomarker values, where higher values indicate a greater likelihood of developing irAEs.
We conclude that all presented models effectively describe the data within the observed range of biomarker values. By 'effectively,' we mean that the optimal model follows an S-shaped curve with an inflection point between the expected centers of the two groups. However, Bayesian models lack predictive power outside the observed range, as meaningful predictions cannot be made without supporting data.
Notably, in Bayesian modeling, the function exhibits a sharp reversal near zero. This fundamental drawback, thoroughly explained in previous sections, arises from the nature of the modeling process and, given the distribution used in our modeling, cannot be mitigated by numerical adjustments. While it could be reduced or even eliminated with alternative distributions, as mentioned at the end of Sec.~\ref{sec:Properties of the Bayesian approach},  this would typically require a greater number of parameters, which, in this case, is neither justified nor feasible with such a limited dataset. 
Therefore, we favor logistic regression, which, despite its simplicity and long-standing use, aligns more closely with the expected function behavior. Furthermore, in the special case where the dispersion of observed data classes is equal, Bayesian modeling reduces to logistic regression.

Beyond selecting the optimal model, an equally critical aspect of our study is understanding the uncertainties associated with these models. While logistic regression demonstrated superior predictive performance within the observed data range, its reliability depends on how well the underlying assumptions hold. To better assess these limitations, we now turn to an in-depth analysis of uncertainty sources and their impact on prediction accuracy.

\paragraph{Quantifying uncertainties:}
The presented models are imperfect mathematical idealizations of reality, containing uncertainties classified as aleatory (data) or epistemic (model). Data uncertainty arises from noise, such as measurement errors, and remains irreducible. 
Model uncertainty stems from limited knowledge, inadequate model structure, or data distribution gaps and can be reduced with more data, complex models, or regularization.
In this study, epistemic uncertainty has the greatest impact, primarily due to data limitations and the model of the data distribution. Model parameters are estimated by fitting to observed data. Key uncertainties affecting our model include:

\begin{enumerate}
    \item Data uncertainty, stemming from variations across scans, noise, and measurement errors.
    \item Model uncertainty, arising from the chosen probabilistic model structure for conditional density functions $p_X(x \mid \ae)$ and $p_X(x \mid \nc)$.
    \item Statistical uncertainty, in estimating parameters of the probabilistic models.
\end{enumerate}

\paragraph{How to use the results:}
Logistic regression is the most widely used modeling approach, as demonstrated by numerous literature reviews on the topic \cite{shipe2019developing, boateng2019review, schober2021logistic}. However, uncertainty analysis, particularly the calculation of confidence intervals or bands, is commonly neglected in the reviewed studies. In this paper, our aim is to address this gap by fitting different binary classification probability models and estimating their uncertainties. The confidence intervals in the resulting model offer valuable insights into the reliability and applicability of the fitted model across different biomarker regimes.

To provide a concrete example, let us consider a hypothetical patient with the value of biomarker $SUV_{\%}$, $x=2$ and discuss the probability of having adverse events based on obtained results.
We analyze the results from different models used to predict the likelihood of adverse events. 
The optimal values represent the best predictions of the model, while the median provides the central value of the predicted distribution, which can be useful for assessing the most likely outcome. Additionally, the confidence interval offers insight into the precision of the predictions and encompasses potential variations of the model, allowing us to evaluate the uncertainty in the predictions. By analyzing these values, we gain a broader understanding of the reliability and applicability of the models, as well as the varying degrees of uncertainty present in certain cases. 

It is important to note that the results are based on $95\%$ confidence intervals, which, in the context of predicting medical outcomes, represent a relatively high standard. This level of uncertainty is often considered a stringent requirement in medical prediction models, as it ensures a higher degree of reliability in the results, but also highlights the inherent variability and complexity involved in predicting clinical events.

 By looking the plots, represented in Result section, concretely at Figures \ref{fig:logit_mvn}, \ref{fig:Gauss_bs} and \ref{fig:Log_normal_bs}, we can collect results for probability of irAE in lungs at biomarker value $x=2$ in Table~\ref{tab:example}.
%
\begin{table}[!htb]

\caption{Collecting results for conditional probability to encounter irAE in lungs, for biomarker value $x=2$. For simpler interpretation probabilities are presented in percentages. $95\%$ confidence intervals $CI = [L_\alpha(x), U_\alpha(x)]$ are formed with significance level $\alpha = 0.05$ as explained by Eq.~\eqref{eq:conf_int}. }

\label{tab:example}
\newcolumntype{C}[1]{>{\centering\arraybackslash}m{#1}}

%\sisetup{round-mode=places,round-precision=3}
\newcommand\perc[2][round-precision = 2]{% default precision: 2
    \SI[round-mode = places,
        scientific-notation = fixed, fixed-exponent = 0,
        output-decimal-marker={.}, #1]{#2e2}{\percent}%
}
\centering
\begin{tabular}{c|C{3.5cm}|C{3.5cm}|C{4cm}|}
Value type &  Logistic regression (\%) & Bayesian-normal model (\%) & Bayesian-lognormal model (\%)\\ \hline
optimal & \perc{0.4233901220524278} & \perc{0.4699474053976193} & \perc{0.4482002564354698} \\
median	& \perc{0.42470951308011934}   & \perc{0.623216859317194}  & \perc{0.5459891837691182} \\
$L_\alpha(x)$   & \perc{0.11738712465726915} & \perc{0.1197384899797259} & \perc{0.18154603727740523} \\
$U_\alpha(x)$	    & \perc{0.8032612583732215}  & \perc{0.9942089585927604} & \perc{0.9247083133161718}
\end{tabular}

\end{table}

The logistic regression model predicts a $42.34\%$ probability of adverse events at the given biomarker value, while the median probability $42.47\%$ is very close to the optimal, indicating that the model's predictions do not vary greatly from the central tendency. The $CI = (11.74\%,80.33\%)$ is very wide, suggesting huge uncertainty in the logistic regression model. 

The Bayesian models, with optimal predictions of $46.99\%$ and $44.82\%$ (for the normal and lognormal distributions, respectively), yield slightly higher probabilities of adverse events compared to the logistic regression model.  Furthermore, the median predictions of $62.32\%$ and $54.60\%$ for the Bayesian-normal and Bayesian-lognormal models are significantly higher than the median of the logistic regression model, reflecting the Bayesian models' tendency to assign greater probabilities to adverse events.

The confidence intervals, $CI = (11.97\%, 99.42\%)$ for the Bayesian-normal model and $CI = (18.15\%, 92.47\%)$ for the Bayesian-lognormal model, are notably wider than those of the logistic regression model. This highlights a higher degree of uncertainty, especially considering the upper bounds reaching almost $100\%$, suggesting that the Bayesian-normal model is more sensitive to variability in the data. While the Bayesian-lognormal model’s predictions are less extreme than the Bayesian-normal model’s, they still exhibit considerable variation when compared to the logistic regression model, which offers more stable predictions.

The broader intervals in the Bayesian models indicate that, while they provide a wider range of potential outcomes, they also capture more variability and uncertainty. 
The logistic regression model provides more conservative and stable predictions with narrower confidence intervals, suitable when lower uncertainty is preferred, but with potentially less sensitivity to higher risk.

These results highlight the trade-off between model complexity, risk sensitivity, and degree of uncertainty, which is crucial when making medical predictions. We can conclude that all models are characterized by a high degree of uncertainty and requiring $95\%$ confidence intervals in such cases may be overly stringent, especially given the limited sample size. This substantial uncertainty is clearly reflected in the presented results, where the wide confidence intervals suggest considerable variability in the predictions. The small sample size further exacerbates this issue, as it limits the model's ability to reliably capture the true underlying distribution and make precise predictions.

While uncertainty quantification provides valuable insights into model reliability, it also highlights the limitations inherent in the current modeling approach. The primary constraints stem from data availability and the choice of biomarkers, both of which influence the robustness of our findings. We highlight two main challenges and limitations encountered during our research in more detail and suggest potential improvements for future research:

\begin{enumerate}

\item The amount of data available when fitting distributions is one of the most significant limiting factors which stems from the specificity of the real-world data we are working with and the inability to expand the modeled dataset significantly. It influences both the choice of model and the optimization of the parameter values within the distribution, leading to the introduction of substantial uncertainties into the modeling process.
A potential avenue for expanding the dataset may allow for the analysis of more complex distributions and the inclusion of additional parameters, leading to a more comprehensive understanding.
We anticipate that the expansion of the database would greatly enhance the model's efficiency, contributing to improved reliability, better fit, and narrower confidence intervals. Currently, the limited dataset poses a significant constraint, making the method somewhat imprecise. 
However, the primary objective of this work is not the results derived from this small dataset but rather the introduction of a novel approach to data access and modeling methodology within the context of cancer treatment practices.

\item In addition, some of the model's imperfections may stem from the choice of biomarker used as the primary predictive feature. Although previous studies \cite{Hribernik2021, Huff2021}, have shown that the biomarker $SUV_\%$ can be used to predict irAEs - specifically, that its increased activity is associated with the occurrence of irAEs, we believe that a significant focus of future research should be on the development of biomarkers. The investigation may explore in more detail the distributions of biomarker values, such as $SUV_\%$ and potential $maxSUV$. Additionally, there is potential for incorporating multiple distinct biomarkers and exploring appropriate aggregation functions to combine them into a unified predictive feature.
We strongly believe that the model can be improved by enhancing its predictive power through the inclusion of more variables in the formation of the biomarker, i.e., the value of the independent variable.

\end{enumerate}



\section{Conclusions}

In this paper, we have implemented two different mathematical approaches for risk quantification of therapeutic process with specific likelihoods, in the treated group of 58 melanoma patients.  We have shown that some of the well-known mathematical methods can be used successfully for risk prediction, with an analysis of the effectiveness and uncertainty of the input data and model. The increased activity of 18F-FDG detected in PET/CT images and measured through a biomarker $SUV_\%$, was shown to be useful in the detection and monitoring of irAE. The mathematical models implemented to quantify information, obtained from 18F-FDG images, should be viewed as powerful tools with great potential to guide personalized treatment strategies in immunotherapy and other cancer therapy treatments.  It could represent the beginning of a new era of data analysis in healthcare practice and should open up opportunities for further work and development in monitoring the therapeutic process, all to avoid potentially life-threatening treatment-related complications.

\section{Acknowledgements}
This research was supported by the Science Fund of the Republic of Serbia, GRANT No 9393, Optimization and prediction in therapy treatments of cancer - OPTIC, by the Ministry of Science, Technological Development and Innovation (Contract No. 451-03-137/2025-03/200156) and the Faculty of Technical Sciences, University of Novi Sad through project “Scientific and Artistic Research Work of Researchers in Teaching and Associate Positions at the Faculty of Technical Sciences, University of Novi Sad 2025” (No. 01-50/295), and by Slovenian Research Agency ARIS, research program P1-0389 and project N1-0197.

\section{Ethical statement}
The study was conducted following the ethical standards defined by the Declaration of Helsinki. It was approved by the Institutional Review Board Committee of both Institutions (Approval number: 2016-0418 in Madison, USA; ERIDKE-0005/2020 in Ljubljana, Slovenia). 

All included subjects consented to being included in the study. At OIL, patients signed informed consent for treatment and consent allowing the usage of their data for scientific purposes. At the UWCCC, the study was approved with a waiver of informed consent.



%\nocite{*}
\section*{References}
\bibliography{refs}

\end{document}