qra/manuscript/optimization.tex

\documentclass{iopjournal}

\input{defs_iop}

\begin{document}

%\articletype{Paper}
\def\journalname{Physics in Medicine \& Biology}
\fancyhead[L]{{\small \sf IOP Publishing}\hspace{5mm} {\it \journalname} {\bf vv} (yyyy) aaaaaa}
\iopheader{\journalname}{dd Month 2026}{dd Month 2026}{Paper}
\fancyhead[R]{Alirezaei {\it et al}\ }

\title{Uncertainty aware risk assessment of metastatic melanoma immunotherapy}

\author{
Zahra Alirezaei$^{1}$\orcid{0000-0002-4094-7042},
Martin Horvat$^1$\orcid{0000-0002-0504-6003},
Marija Delić$^2$\orcid{0000-0002-3335-7546},
Katja Strašek$^1$\orcid{0009-0004-6115-1512},
Daniel Huff$^3$\orcid{0000-0001-9792-4119} and
Robert Jeraj$^{1,3}$\orcid{0000-0002-2192-2931}
}


\affil{$^1$Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana, Slovenia}

\affil{$^2$Faculty of Technical Sciences, University of Novi Sad, Serbia}

\affil{$^3$Department of Medical Physics, School of Medicine and Public Health, University of Wisconsin--Madison, Madison, WI, USA}

\affil{$^*$Author to whom any correspondence should be addressed.}

\email{zahra.alirezaei@fmf.uni-lj.si}

\keywords{immune checkpoint inhibitors, pneumonitis, PET/CT, quantitative imaging biomarkers, uncertainty quantification, logistic regression, Bayesian modeling, metastatic melanoma}

\begin{abstract}

\textit{Objective}: Immune checkpoint inhibitors (ICIs) have improved outcomes in metastatic melanoma but can also cause immune-related adverse events, including pneumonitis. As quantitative imaging biomarkers move from identifying associations with toxicity toward modeling patient risk, an important question arises: how certain are the estimated risks when the available data are limited, and the biomarker itself is subject to measurement variability? This question is particularly relevant in small and imbalanced datasets, where modeling assumptions and individual observations may also affect the estimated risk. This study, therefore, aims to develop an uncertainty-aware framework to quantitatively assess pneumonitis risk using $^{18}$F-FDG PET/CT-derived lung biomarkers.

\textit{Approach}: For each patient, the imaging biomarker was defined as the maximum longitudinal lung  \(\SUV_{95\%}\) extracted from serial $^{18}$F-FDG PET/CT examinations. Pneumonitis risk was modeled using standard logistic regression, a monotone cubic logistic extension, and a constrained Bayesian generative model incorporating monotonicity and tail constraints. Uncertainty was assessed for model parameters, fitted risk functions, and derived quantities using complementary uncertainty quantification approaches, including Wald-based intervals and nonparametric and parametric bootstrap procedures with 1000 replicates. Measurement-noise propagation was used to examine how uncertainty in the imaging biomarker affected the estimated risk, while sensitivity analysis examined the effects of local changes in model parameters. Analyses of the full dataset and a trimmed dataset excluding an identified high-leverage observation from the normal-control group were used to assess the influence of this observation on model estimates and uncertainty.

\textit{Main results}: All three model formulations produced a monotone-increasing relationship between the lung biomarker and estimated pneumonitis risk. The analyses showed that uncertainty in model parameters, fitted risk functions, and derived quantities did not necessarily follow the same pattern. Comparison of the full and trimmed datasets revealed differences in the fitted risk functions and their uncertainty after removing the identified observation. Measurement-noise propagation showed how uncertainty in the imaging biomarker propagated to the estimated risk, while sensitivity analysis showed differences in the response of the derived quantities to changes in individual model parameters. The constrained Bayesian model produced a smooth monotone risk function while incorporating explicit constraints on monotonicity and tail behavior.

\textit{Significance}: Moving from estimating a biomarker--risk relationship to uncertainty-aware risk modeling requires consideration of uncertainty beyond model fit alone. The proposed framework combines model estimation, uncertainty quantification, measurement-noise propagation, and sensitivity analysis to examine uncertainty in the resulting risk estimates. This provides a methodological basis for assessing $^{18}$F-FDG PET/CT-based risk models under limited and imbalanced data conditions. External validation and evaluation of clinical utility are required before clinical application.

\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

Immune checkpoint inhibitors (ICIs) have substantially improved outcomes in metastatic melanoma and other malignancies \citep{wolchok2017overall, wolchok2022long}. However, they can also induce immune-related adverse events (irAEs), which range from mild to life-threatening \citep{wang2018fatal, puzanov2017managing}. Among these, immune checkpoint inhibitor-related pneumonitis is of particular clinical concern because of its potential severity and its impact on treatment continuation. Identifying imaging-derived biomarkers associated with pneumonitis risk may therefore contribute to the quantitative assessment of toxicity risk during ICI therapy.

Whole-body $^{18}$F-FDG PET/CT ($^{18}$F-fluorodeoxyglucose positron emission tomography/computed tomography) is routinely used for disease assessment and monitoring treatment response during ICI therapy \citep{filippi2022_diagnostics, aide2022_cancerimaging, gandy2020_bjr}. Beyond tumor metabolism, $^{18}$F-FDG uptake also reflects inflammatory processes, making PET/CT a potential source of quantitative imaging biomarkers (QIBs) related to immune-mediated toxicity \citep{iravani2020_ejnmmi, nadaraja2024_cti}. Previous studies have reported associations between PET-derived biomarkers and immune-related adverse events, including lung uptake measures associated with pneumonitis \citep{Huff2021, Hribernik2022}. These findings motivate the development of quantitative models that relate imaging-derived biomarker values to estimated adverse-event risk.

Despite these advances, translating PET-derived biomarkers into quantitative risk estimates remains challenging. Estimating a biomarker--risk relationship is only one component of predictive risk modeling; the uncertainty associated with that relationship must also be considered. This is particularly important in small and imbalanced clinical datasets, where apparently well-defined risk estimates may coexist with uncertainty arising from parameter estimation, limited sampling, individual observations, or uncertainty in the measured biomarker. Uncertainty quantification has been recognized as an important component of predictive modeling in medical applications, particularly when model outputs may ultimately inform clinical decision-making \citep{Begoli2019, Tomasev2021}. Uncertainty is also an important consideration in medical imaging, where both the imaging process and subsequent model-based analysis can contribute to uncertainty in the resulting estimates \citep{Lambert2024}. Consequently, the assessment of an imaging-biomarker risk model should consider both the uncertainty associated with the fitted model and the propagation of biomarker uncertainty into the resulting risk estimates.

Uncertainty in a fitted risk model can be characterized at different levels, which do not necessarily provide equivalent information. Parameter intervals describe the precision of model-specific parameter estimates, whereas confidence bands for the fitted risk function characterize the precision of estimated risk across the biomarker range. Uncertainty can also be evaluated for quantities derived from the fitted risk function. Moreover, greater model flexibility may improve representation of the observed data without necessarily improving parameter identifiability or reducing uncertainty in the resulting risk estimates. An uncertainty-aware assessment should therefore consider these complementary aspects of model behavior rather than relying on a single measure of model fit or uncertainty.

This study aims to develop and evaluate an uncertainty-aware quantitative risk assessment (QRA) framework for estimating the probability of immune checkpoint inhibitor-related pneumonitis as a function of $^{18}$F-FDG PET/CT-derived lung biomarkers. We compare alternative probabilistic model formulations, including standard logistic regression, a monotone cubic logistic extension, and a constrained Bayesian generative model, and examine their behavior on both the full and trimmed datasets obtained by excluding the identified normal-control observation. Uncertainty is assessed for model parameters, fitted risk functions, and derived quantities using complementary uncertainty quantification approaches. In addition, we investigate how measurement noise in the imaging biomarker propagates through the fitted risk models and how local changes in model parameters affect the derived quantities through sensitivity analysis.

The proposed framework brings these analyses together to examine the fitted risk models under the limited and imbalanced data conditions considered in this study. The present work focuses on methodological development and uncertainty-aware risk assessment; external validation and evaluation of clinical utility are required before clinical application.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Methods}
\label{sec:methods}

\subsection{Patient data and biomarker extraction}
\label{sec:data_biomarker}

We retrospectively analyzed $^{18}$F-FDG PET/CT imaging and clinical data from 58 patients with metastatic melanoma treated with immune checkpoint inhibitors (anti-CTLA-4 and/or anti-PD-1) at the Institute of Oncology Ljubljana (Slovenia) and the University of Wisconsin Carbone Cancer Center (Madison, USA). The adverse event ($\AEgrp$) of interest was immune checkpoint inhibitor-related pneumonitis, ascertained by chart review. Patients were assigned to the normal-control group ($\NCgrp$) if pneumonitis did not occur, and to the AE group ($\AEgrp$) if it did \citep{Hribernik2022}.

Image acquisition, reconstruction, and lung parenchyma segmentation followed the previously published protocol \citep{Hribernik2022}. Lung segmentation was performed using the DeepMedic convolutional neural network architecture \citep{Kamnitsas2017DeepMedic}. Quantitative imaging biomarkers were defined as organ-wise $^{18}$F-FDG PET/CT uptake percentiles, denoted $\SUV_{p\%}$, computed within the segmented lung volume \citep{Huff2021}. Consistent with prior analyses, we selected the upper-tail percentile $\SUV_{95\%}$ as the primary biomarker.

For each patient \(i\), we used the maximum lung \(\SUV_{95\%}\) across eligible pre-event or early-treatment visits as the patient-level biomarker, following the biomarker aggregation strategy described previously \citep{Huff2021}.
% -------------------------------------------------------------------------
\subsection{Modeling framework and data variants}
\label{sec:framework}

We consider a random pair \((X,Y)\), where \(X \ge 0\) denotes the patient-level \(^{18}\mathrm{F}\)-FDG PET/CT-derived lung biomarker \(\SUV_{95\%}\), as defined in the previous section, and \(Y \in \{\NCgrp,\AEgrp\}\) denotes pneumonitis status.

Our goal is to model the individual risk of an adverse event (\(\AEgrp\)) as a function of the observed biomarker value. In this setting, the risk at the biomarker level \(x\) is defined as
%
\begin{equation}
    \label{eq:target_main}
    \Prob(Y=\AEgrp \mid X=x),
\end{equation}
%
that is, the conditional probability that a patient belongs to the adverse-event group given \(X=x\). We represent this risk by a parametric function:
%
\begin{equation}
    \label{eq:riskcurve_main}
    p(x;\bTheta) \equiv \Prob_{\bTheta}(Y=\AEgrp \mid X=x),
\end{equation}
%
where \(\bTheta\) denotes the model parameter vector.

The model formulations considered in this study were evaluated with respect to the following properties of the resulting risk function:
%
\begin{enumerate}
    \item \(p(x;\bTheta)\in[0,1] \quad \forall x \ge 0\), so that the predicted risk remains within the valid probability range;
    \item \(p'(x) \ge 0\), corresponding to a monotonically non-decreasing relationship between the biomarker and predicted risk;
    \item
    \[
    \lim_{x \to 0^+} p(x;\bTheta) = 0,
    \qquad
    \lim_{x \to \infty} p(x;\bTheta) = 1,
    \]
    corresponding to lower predicted risk at the lower biomarker boundary and increasing predicted risk toward unity as the biomarker becomes large.
\end{enumerate}
%
These properties provide a common basis for examining the behavior of the alternative risk models. Their implementation, however, depends on the model formulation: probability bounds follow from the probabilistic model structure, whereas monotonicity and tail behavior are imposed explicitly where required, as described in the corresponding model sections below.

Each patient contributes one pair \((x_i,y_i)\), where \(x_i\) is the patient-level biomarker value and \(y_i\) is the associated group label. For \(n\) patients, the resulting dataset is
%
\begin{equation}
    \label{eq:dataset_main}
    \mathcal{D} = \{(x_i, y_i)\}_{i=1}^{n}.
\end{equation}
%

Analyses were performed on two related datasets summarized in Table~\ref{tab:datasets}. The primary FULL dataset contained all 58 patients. A TRIM dataset was additionally constructed by excluding one \(\NCgrp\) observation with lung \(\SUV_{95\%}=2.48\), which lies in the upper tail of the \(\NCgrp\) biomarker distribution and overlaps with the \(\AEgrp\) range. Given the small and imbalanced dataset, the TRIM analysis was used to examine the sensitivity of the fitted risk models and their uncertainty to the removal of this observation. The FULL dataset remained the primary dataset for interpretation.

\begin{table}[!htbp]
\centering
\caption{Composition of the analyzed metastatic melanoma datasets used for PET-derived lung \(\SUV_{95\%}\) biomarker modeling. AE denotes patients who developed immune checkpoint inhibitor-related pneumonitis, and NC denotes patients without pneumonitis.}
\label{tab:datasets}
\begin{tabular}{lccc}
\hline
Dataset & Total & AE & NC \\
\hline
FULL & 58 & 5 & 53 \\
TRIM & 57 & 5 & 52 \\
\hline
\end{tabular}
\end{table}
% -------------------------------------------------------------------------
\subsection{Risk modeling approaches}
\label{sec:risk_models}

\subsubsection{Logistic regression models}
\label{sec:logistic}

As a baseline discriminative approach, we model the risk of AE using logistic regression (LR) \citep{McCullaghNelder1989, Hastie2009Elements}. Here, risk refers to the conditional probability of an adverse event given the biomarker value \(X=x\). The logistic regression model is given by
%
\begin{equation}
  \label{eq:logistic_main}
  p(x;\bTheta)
  =
  \expit\left(\eta_{\mathrm{LR}}(g(x);\bTheta)\right),
  \qquad
  \expit(u)=\frac{1}{1+e^{-u}},
\end{equation}
%
where the log-odds are modeled as a linear function of the transformed predictor,
%
\begin{equation}
  \eta_{\mathrm{LR}}(t;\bTheta)
  =
  \theta_0+\theta_1 t,
\end{equation}
%
with two predictor parameterizations:
%
\begin{equation}
  g(x)=x \quad \text{(RAW)},
  \qquad
  g(x)=\log x \quad \text{(LOG)}.
\end{equation}

\paragraph{Monotone cubic logistic extension.}

As a more flexible shape-constrained comparator, we also considered a monotone cubic logistic model of the form
%
\begin{equation}
  p(x;\bTheta)
  =
  \expit\left(\eta_{\mathrm{cubic\,LR}}(g(x);\bTheta)\right),
\end{equation}
%
with cubic log-odds
%
\begin{equation}
  \eta_{\mathrm{cubic\,LR}}(t;\bTheta)
  =
  \theta_0+\theta_1 t+\theta_2 t^2+\theta_3 t^3.
\end{equation}
%
Monotonicity was enforced through a constrained reparameterization that ensures
\(\eta_{\mathrm{cubic\,LR}}'(t;\bTheta)\ge 0\) over the specified predictor domain. This allows greater functional flexibility while maintaining a non-decreasing risk function. Full technical details of the parameterization and monotonicity constraint are provided in Appendix~\ref{sec:app_logistic}.

\paragraph{Model fitting.}

For the standard logistic models, the parameters
\(\bTheta=(\theta_0,\theta_1)\) were estimated by maximum likelihood. Assuming independent Bernoulli observations, the log-likelihood is
%
\begin{equation}
\ell_{\mathrm{LR}}(\bTheta\mid\mathcal D)
=
\sum_{i=1}^{n}
\left[
y_i\log p(x_i;\bTheta)
+
(1-y_i)\log\left(1-p(x_i;\bTheta)\right)
\right],
\label{eq:bernoulli_llf_main}
\end{equation}
%
where \(p(x_i;\bTheta)\) denotes the predicted probability of an adverse event for patient \(i\). The maximum likelihood estimator (MLE) is
%
\begin{equation}
\hat{\bTheta}
=
\argmax_{\bTheta}
\ell_{\mathrm{LR}}(\bTheta\mid\mathcal D).
\label{eq:mle_main}
\end{equation}
%
For the standard logistic models, optimization was performed using the L-BFGS algorithm \citep{Nocedal2006}. For the monotone cubic logistic model, parameter estimation used a two-stage strategy consisting of differential evolution for global exploration followed by L-BFGS refinement \citep{Storn1997, Virtanen2020, Nocedal2006}. Under standard regularity conditions, the MLE is asymptotically normally distributed, providing the basis for the asymptotic uncertainty quantification described below. Additional mathematical details, including likelihood derivatives and covariance expressions, are provided in Appendix~\ref{sec:app_logistic}.
% -------------------------------------------------------------------------
\subsubsection{Bayesian probability modeling}
\label{sec:bayesian_model}

Bayes' theorem provides a framework for modeling risk by combining observed data with assumptions about the underlying distributions \citep{Bishop2006, Gelman2013}. In this context, risk is defined as the conditional probability of an adverse event (AE) given the observed biomarker value \(X\). Let \(Y\in\{\AEgrp,\NCgrp\}\) denote whether a patient belongs to the AE or normal-control (NC) group. Using Bayes' theorem,
%
\begin{equation}
P(Y=\AEgrp \mid X=x)
=
\frac{
p_{X\mid Y}(x\mid\AEgrp)\,P_Y(\AEgrp)
}{
p_X(x)
},
\end{equation}
%
where \(p_{X\mid Y}(x\mid y)\) is the group-conditional biomarker density, \(P_Y(y)\) is the marginal probability of group \(y\), and \(p_X(x)\) is the marginal biomarker density. By the law of total probability,
%
\begin{equation}
p_X(x)
=
p_{X\mid Y}(x\mid\AEgrp)P_Y(\AEgrp)
+
p_{X\mid Y}(x\mid\NCgrp)P_Y(\NCgrp).
\end{equation}
%
Substitution gives
%
\begin{equation}
P(Y=\AEgrp \mid X=x)
=
\frac{
p_{X\mid Y}(x\mid\AEgrp)\,P_Y(\AEgrp)
}{
p_{X\mid Y}(x\mid\AEgrp)\,P_Y(\AEgrp)
+
p_{X\mid Y}(x\mid\NCgrp)\,P_Y(\NCgrp)
}.
\label{eq:bayes_cond_prob}
\end{equation}
%
Thus, the risk \(P(Y=\AEgrp \mid X=x)\) depends on the AE prevalence and the two group-conditional densities evaluated at \(x\). Since
%
\[
P_Y(\NCgrp)+P_Y(\AEgrp)=1,
\]
%
the risk is determined by three components: \(P_Y(\AEgrp)\), \(p_{X\mid Y}(x\mid\AEgrp)\), and \(p_{X\mid Y}(x\mid\NCgrp)\). Consequently, the choice of distribution families for the two group-conditional densities affects not only within-group fit but also the shape, monotonicity, and asymptotic behavior of the resulting risk function. Further details on the distributional choices and their implications are provided in Appendix~\ref{sec:app_families}.

\paragraph{Distributional model selection.}

To model the risk \(P(Y=\AEgrp \mid X=x)\), we adopt a generative formulation in which the group-conditional densities and prevalence are parameterized explicitly. Several candidate positive-support distribution families were evaluated before selecting the final specification.

Distribution selection was guided not only by empirical within-group fit but also by the behavior of the induced risk function. In particular, some combinations of class-conditional distributions, including Normal and log-normal models, can produce non-monotone or unstable risk behavior even when the individual group fits appear adequate. Representative examples and mathematical details are provided in Appendix~\ref{sec:app_families}. The final distributional model was therefore selected by considering within-group goodness of fit, stability, monotonicity of the induced risk function, and its asymptotic tail behavior.

Based on these considerations, the final model uses a Beta-prime distribution \citep{Johnson1995} for the AE group and a Gamma distribution \citep{Johnson1994} for the NC group. Explicit probability density functions and the alternative candidate families considered are provided in Appendix~\ref{sec:app_families_candidates}.
%
\begin{align}
    X \mid (Y=\AEgrp) &\sim \mathrm{BetaPrime}(a,b,s), \\
    X \mid (Y=\NCgrp) &\sim \mathrm{Gamma}(k,\vartheta).
\end{align}
%
This specification accommodates positive, right-skewed biomarker values and, together with the parameter constraints introduced below, yields a monotone risk function with the specified asymptotic behavior. The model is therefore structured so that monotonicity and tail behavior of the resulting risk function can be imposed through analytic parameter constraints.

\paragraph{Model assumptions and constraints.}

The constrained Bayesian model is constructed to satisfy three properties: valid probabilities, a monotonically non-decreasing risk function, and the specified tail behavior.

For the selected Beta-prime--Gamma distributional model, the asymptotic behavior of the log-odds implies that \(a>k\) is required for the risk to approach zero as \(x\to0^+\). The Gamma scale parameter satisfies \(\vartheta>0\), and the exponential decay of the Gamma density relative to the Beta-prime upper tail gives the corresponding upper-tail behavior as \(x\to\infty\). These conditions alone, however, do not guarantee global monotonicity. To enforce monotonicity, the derivative of the log-odds is constrained to be non-negative for all \(x>0\). This leads to the analytically derived admissible bound
%
\begin{equation}
  \vartheta
  \le
  \frac{s}
  {\left(\sqrt{a+b}-\sqrt{a-k}\right)^2},
  \qquad a>k,
\end{equation}
%
which defines the feasible parameter region for which the resulting risk function is monotone. The derivation of this condition is provided in Appendix~\ref{sec:app_bayes_conditions}.

\paragraph{Model fitting.}

Using Eq.~\eqref{eq:bayes_cond_prob}, we define the risk function as
%
\begin{equation}
p(x;\bTheta)
\equiv
P(Y=\AEgrp \mid X=x),
\end{equation}
%
with model parameter vector
%
\[
\bTheta=(\omega,a,b,s,k,\vartheta),
\]
%
where \((a,b,s)\) are the parameters of the Beta-prime distribution for the \(\AEgrp\) group, \((k,\vartheta)\) are the parameters of the Gamma distribution for the \(\NCgrp\) group, and
%
\[
\omega=P_Y(\AEgrp)
\]
%
is the AE prevalence. The conditional probability for either group can then be written as
%
\begin{equation}
p(y\mid x,\bTheta)
=
\begin{cases}
1-p(x;\bTheta), & y=\NCgrp,\\
p(x;\bTheta),   & y=\AEgrp.
\end{cases}
\end{equation}
%
The log-posterior objective for the constrained Bayesian model is
%
\begin{equation}
L_{\mathrm{CB}}(\bTheta\mid\mathcal{D})
=
\sum_{i=1}^{n}
\log p(y_i\mid x_i,\bTheta)
+
\log p(\bTheta),
\label{eq:bayes_objfun_main}
\end{equation}
%
where \(\mathcal D\) denotes the observed dataset introduced in Section~\ref{sec:framework}, and \(p(\bTheta)\) denotes the prior density on the model parameters. A detailed breakdown of the objective function is provided in Appendix~\ref{sec:app_bayes_objfun}.

The first term in Eq.~\eqref{eq:bayes_objfun_main} is the conditional log-likelihood of the observed group labels given the biomarker values, whereas the second term represents the contribution of the parameter priors. The parameter estimate is obtained by maximizing the log-posterior over the feasible parameter region:
%
\begin{equation}
\hat{\bTheta}_{\mathrm{MAP}}
=
\argmax_{\bTheta\in\Omega}
L_{\mathrm{CB}}(\bTheta\mid\mathcal D),
\end{equation}
%
where \(\Omega\) denotes the feasible parameter region defined by the structural constraints. The resulting estimate is the maximum a posteriori (MAP) estimate.

Rather than optimizing the constrained parameters \(\bTheta\) directly, optimization was performed in an unconstrained internal parameter space. Parameters were mapped to \(\bTheta\) using smooth transformations: the prevalence \(\omega\) was mapped to \((0,1)\) using a sigmoid transformation, positive parameters were obtained using softplus transformations, and the Gamma scale parameter \(\vartheta\) was expressed relative to its analytically derived admissible upper bound. A formal description of these transformations is provided in Appendix~\ref{sec:app_bayes_transf}.

This reparameterization ensures that evaluated parameter vectors satisfy the structural constraints by construction, converting the constrained optimization problem into an unconstrained one. The MAP formulation combines the conditional Bernoulli likelihood with weakly informative parameter priors.

\subsection{Model evaluation}
\label{sec:model_evaluation}

Model evaluation was performed using the Bernoulli log-likelihood (LLF), the Akaike information criterion (AIC) \citep{Akaike1974}, and the Bayesian information criterion (BIC) \citep{Schwarz1978}. For a fitted model with log-likelihood \(\ell(\hat{\bTheta})\), \(q\) estimated parameters, and \(n\) observations, the information criteria are defined as
%
\begin{align}
\mathrm{AIC} &= 2q - 2\,\ell(\hat{\bTheta}),\\
\mathrm{BIC} &= q\log n - 2\,\ell(\hat{\bTheta}).
\end{align}
%
For the standard logistic regression models, \(\ell(\hat{\bTheta})\) corresponds to the Bernoulli log-likelihood evaluated at the maximum likelihood estimate. For the constrained Bayesian model, the Bernoulli log-likelihood was evaluated at the MAP estimate to compare with the standard logistic regression models; the prior contribution was not included in the LLF used to calculate AIC and BIC. The parameter count was \(q=2\) for each standard logistic regression model and \(q=6\) for the constrained Bayesian model. AIC and BIC were not reported for the monotone cubic logistic model.
% -------------------------------------------------------------------------
\subsection{Risk function characteristics}
\label{sec:risk_characteristics_method}

To summarize and compare the fitted risk functions
\(p(x;\hat{\bTheta})\) across model formulations, we considered two derived quantities:

\begin{itemize}

\item Midpoint \(x_{50}\) -- the biomarker value at which the predicted risk reaches 50\%:

\begin{equation}
p(x_{50};\hat{\bTheta})=\frac{1}{2}.
\label{eq:x50}
\end{equation}

\item Midpoint slope \(s_{50}\) -- the slope of the fitted risk function evaluated at \(x_{50}\):

\begin{equation}
s_{50}
=
\left.
\frac{\partial}{\partial x}
p(x;\hat{\bTheta})
\right|_{x=x_{50}}.
\label{eq:s50}
\end{equation}

\end{itemize}

The midpoint \(x_{50}\) provides the biomarker value corresponding to a predicted risk of 50\% and is analogous to the median-effect point (ED$_{50}$) commonly used in dose--response modeling \citep{Ritz2015DoseResponse}. The midpoint slope \(s_{50}\) provides the corresponding rate of change in predicted risk with respect to the biomarker at \(x_{50}\). Together, these quantities provide common summaries of the fitted risk function that can be evaluated across model formulations without requiring direct comparison of model-specific parameters. Uncertainty in \(x_{50}\) and \(s_{50}\) was therefore evaluated alongside uncertainty in the model parameters and fitted risk functions.

% -----------------------------------------------------------------------
\subsection{Sensitivity Analysis}
\label{sec:sensitivity_elasticity}

To assess how model parameters influence the derived risk-function characteristics, we evaluated local sensitivities and relative sensitivities (elasticities) \citep{Saltelli2008} for the midpoint \(x_{50}\) and the corresponding midpoint slope \(s_{50}\).

Let \(\bTheta\) denote the vector of model parameters. For the standard logistic models,
\(\bTheta=(\theta_0,\theta_1)\), whereas for the constrained Bayesian model,
\(\bTheta=(\omega,a,b,s,k,\vartheta)\).
In the following, \(\theta_j\) denotes a generic component of \(\bTheta\), irrespective of its specific name in a given model.

For a risk-function characteristic \(q\), where \(q\in\{x_{50},\,s_{50}\}\), the local sensitivity with respect to parameter \(\theta_j\) is defined as
%
\begin{equation}
\frac{\partial q}{\partial \theta_j}.
\label{eq:local_sensitivity}
\end{equation}
%
This quantity measures the rate of change in a risk-function characteristic in response to a small change in an individual model parameter.

Elasticity provides a dimensionless measure of relative parameter influence and is defined as
%
\begin{equation}
\mathcal{E}_{\theta_j}(q)
=
\frac{\partial q}{\partial \theta_j}
\frac{\theta_j}{q}.
\label{eq:elasticity}
\end{equation}
%
Elasticities facilitate comparison of parameter effects across different scales and units by expressing the change in a risk-function characteristic relative to proportional changes in the corresponding model parameter \citep{Saltelli2008}.

When analytical derivatives were unavailable, local sensitivities and elasticities were approximated using central finite differences \citep{Nocedal2006}. Sensitivities were evaluated with respect to all parameters in \(\bTheta\). For the constrained Bayesian model, parameterizations were transformed where necessary to improve numerical stability during numerical differentiation.
% -------------------------------------------------------------------------
\subsection{Uncertainty quantification}
\label{sec:uq}

Uncertainty quantification (UQ) was used to assess uncertainty in the estimated risk functions and their derived characteristics. Two sources of uncertainty were considered: \emph{epistemic uncertainty}, associated with estimation of the statistical model from a finite dataset, and \emph{aleatoric uncertainty}, arising from variability in the quantitative imaging biomarker. The latter represents the component of aleatoric uncertainty considered in this study. Epistemic uncertainty was assessed using confidence intervals for the model parameters, confidence bands for the fitted risk functions, and confidence intervals for derived characteristics of the risk functions. Aleatoric uncertainty was assessed separately by propagating biomarker perturbations through the fitted models. The methods used for these two components of the UQ analysis are described below.
% ----------------------------------
\subsubsection{Confidence intervals and confidence bands}
\label{sec:uq:intervals}

Uncertainty associated with model estimation was quantified using four complementary approaches: analytical approximation, Monte Carlo approximation, nonparametric bootstrap (NPBS), and parametric bootstrap (PBS). The analytical approximation, NPBS, and PBS were applied to both the logistic regression and constrained Bayesian models, whereas the Monte Carlo approximation was applied only to logistic regression. Uncertainty was evaluated for the model parameters, the fitted risk function \(p(x;\boldsymbol{\Theta})\), and the derived risk-function characteristics applicable to each model. Here, \(x_{50}\) denotes the biomarker value at which the predicted risk equals 0.5, and \(s_{50}\), where applicable, denotes the local slope of the fitted risk function at \(x_{50}\).

\paragraph{Wald approximation}

The Wald approximation provides a local characterization of parameter uncertainty based on the curvature of the fitted objective function \citep{CasellaBerger2002}. Near an optimum, a sufficiently smooth objective function can be approximated by a quadratic function of the model parameters. This local quadratic approximation leads to a Gaussian representation of uncertainty around the fitted parameter vector, with the covariance determined by the inverse curvature at the optimum. A sharply curved objective therefore corresponds to smaller local parameter uncertainty, whereas weak curvature indicates greater uncertainty in the corresponding parameter directions.

The same principle was applied to logistic regression (LR) and the constrained Bayesian (CB) model. The two models differ in the objective function used for parameter estimation: LR parameters are obtained by maximizing the Bernoulli log-likelihood, whereas CB parameters are obtained by maximizing the log-posterior. Their respective estimates are

\begin{align}
\hat{\boldsymbol{\Theta}}_{\mathrm{LR}}
&=
\arg\max_{\boldsymbol{\Theta}_{\mathrm{LR}}}
\ell_{\mathrm{LR}}
\left(
\boldsymbol{\Theta}_{\mathrm{LR}}
\mid \mathcal{D}
\right),
\label{eq:lr_mle}
\\
\hat{\boldsymbol{\Theta}}_{\mathrm{CB}}
&=
\arg\max_{\boldsymbol{\Theta}_{\mathrm{CB}}}
\log p
\left(
\boldsymbol{\Theta}_{\mathrm{CB}}
\mid \mathcal{D}
\right),
\label{eq:cb_map}
\end{align}

where $\ell_{\mathrm{LR}}$ denotes the Bernoulli log-likelihood and
$p(\boldsymbol{\Theta}_{\mathrm{CB}}\mid\mathcal{D})$ denotes the posterior distribution. Thus, although the objective functions differ, the local uncertainty calculation follows the same curvature-based principle.

For either model, let $\boldsymbol{\psi}_{m}$ denote the parameter vector in the space in which the local curvature is evaluated, with

\begin{equation}
\boldsymbol{\psi}_{m}
=
\begin{cases}
\boldsymbol{\Theta}_{\mathrm{LR}}, & m=\mathrm{LR},\\
\boldsymbol{\phi}, & m=\mathrm{CB},
\end{cases}
\label{eq:wald_parameter_space}
\end{equation}

where $\boldsymbol{\phi}$ denotes the unconstrained parameterization of the CB model. The local covariance matrix is approximated by

\begin{equation}
\hat{\boldsymbol{\Sigma}}_{\psi_m}
=
\mathbf{H}_{m}^{-1}
\left(
\hat{\boldsymbol{\psi}}_{m}
\right),
\qquad
m\in\{\mathrm{LR},\mathrm{CB}\},
\label{eq:wald_cov_general}
\end{equation}

where $\mathbf{H}_{\mathrm{LR}}$ is the observed Hessian of the negative log-likelihood and $\mathbf{H}_{\mathrm{CB}}$ is the observed Hessian of the negative log-posterior, each evaluated at the corresponding fitted parameter vector. Accordingly, the local Gaussian approximation can be expressed in the common form

\begin{equation}
\boldsymbol{\psi}_{m}
\mathrel{\dot{\sim}}
\mathcal{N}
\left(
\hat{\boldsymbol{\psi}}_{m},
\hat{\boldsymbol{\Sigma}}_{\psi_m}
\right),
\qquad
m\in\{\mathrm{LR},\mathrm{CB}\}.
\label{eq:wald_general}
\end{equation}

For LR, $\boldsymbol{\psi}_{\mathrm{LR}}=\boldsymbol{\Theta}_{\mathrm{LR}}$ directly. For CB, uncertainty represented in the unconstrained space is mapped to the original parameterization through
$\boldsymbol{\Theta}_{\mathrm{CB}}=T(\boldsymbol{\phi})$.

This local Gaussian representation provides the common basis for uncertainty propagation in both model formulations. In the Wald approach, covariance is propagated to the fitted risk function $p(x;\boldsymbol{\Theta}_{m})$ and to the derived characteristics $x_{50}$ and $s_{50}$ using the first-order delta method \citep{CasellaBerger2002}. The resulting intervals therefore represent a first-order approximation to the uncertainty in the model parameters, fitted risk estimates, and derived characteristics. Model-specific mathematical derivations are provided in Appendices~\ref{sec:app_ci_logistic} and~\ref{sec:app_ci_bayes}.

\paragraph{Monte Carlo approximation (MCA)}

The Monte Carlo approximation (MCA) provides a numerical alternative to the first-order uncertainty propagation used in the Wald approximation. Both approaches start from the same local Gaussian representation of parameter uncertainty obtained from the fitted model and its Hessian. However, whereas the Wald approach propagates the parameter covariance analytically using a first-order Taylor expansion, MCA repeatedly samples parameter vectors from the corresponding Gaussian approximation and evaluates the complete nonlinear risk model for each realization \citep{CasellaBerger2002, RobertCasella2004}. This allows uncertainty to be propagated to the fitted risk function and its derived characteristics without linearizing their dependence on the model parameters.

For logistic regression (LR), Monte Carlo realizations are sampled directly in the model parameter space,

\begin{equation}
\boldsymbol{\Theta}_{\mathrm{LR}}^{(j)}
\sim
\mathcal{N}
\left(
\hat{\boldsymbol{\Theta}}_{\mathrm{LR}},
\hat{\boldsymbol{\Sigma}}_{\mathrm{LR}}
\right),
\qquad
j=1,\ldots,M.
\label{eq:mc_sampling_lr}
\end{equation}

For the constrained Bayesian (CB) model, direct Gaussian sampling in the original parameter space may generate parameter combinations that violate the model constraints. The local Gaussian approximation is therefore constructed in the unconstrained parameter space. Let $\boldsymbol{\phi}$ denote the unconstrained parameter vector and let
$T(\boldsymbol{\phi})$ denote the transformation from the unconstrained to the original CB parameter space. Monte Carlo realizations are generated according to

\begin{align}
\boldsymbol{\phi}^{(j)}
&\sim
\mathcal{N}
\left(
\hat{\boldsymbol{\phi}},
\hat{\boldsymbol{\Sigma}}_{\phi}
\right),
\label{eq:mc_sampling_cb_phi}
\\
\boldsymbol{\Theta}_{\mathrm{CB}}^{(j)}
&=
T\left(\boldsymbol{\phi}^{(j)}\right),
\qquad
j=1,\ldots,M,
\label{eq:mc_transform_cb}
\end{align}

where $\hat{\boldsymbol{\phi}}$ is the fitted CB parameter vector expressed in the unconstrained space and $\hat{\boldsymbol{\Sigma}}_{\phi}$ is obtained from the inverse Hessian of the negative log-posterior with respect to $\boldsymbol{\phi}$. This construction preserves the parameter constraints when the Monte Carlo realizations are mapped back to the original CB parameterization.

For both models, the resulting parameter realizations are propagated through their respective risk functions,

\begin{equation}
p_{m}^{(j)}(x)
=
p
\left(
x;
\boldsymbol{\Theta}_{m}^{(j)}
\right),
\qquad
m\in\{\mathrm{LR},\mathrm{CB}\},
\qquad
j=1,\ldots,M.
\label{eq:mc_curve}
\end{equation}

A total of $M=10\,000$ realizations are generated for each model. The resulting ensemble of risk functions is used to characterize uncertainty across the biomarker range, with pointwise 95\% uncertainty bands obtained from the empirical 2.5th and 97.5th percentiles of
$p_{m}^{(j)}(x)$ at each biomarker value.

The same Monte Carlo realizations are used to propagate parameter uncertainty to the derived risk-function characteristics. For each realization, $x_{50,m}^{(j)}$ is obtained from

\begin{equation}
p
\left(
x_{50,m}^{(j)};
\boldsymbol{\Theta}_{m}^{(j)}
\right)
=
0.5,
\label{eq:mc_x50}
\end{equation}

and the corresponding midpoint slope is evaluated as

\begin{equation}
s_{50,m}^{(j)}
=
\left.
\frac{\partial
p\left(x;\boldsymbol{\Theta}_{m}^{(j)}\right)}
{\partial x}
\right|_{x=x_{50,m}^{(j)}}.
\label{eq:mc_s50}
\end{equation}

The 95\% intervals for $x_{50}$ and $s_{50}$ are obtained from the empirical 2.5th and 97.5th percentiles of their respective Monte Carlo distributions. In contrast to first-order analytical propagation, MCA evaluates the nonlinear mapping from parameter uncertainty to the quantities of interest directly for every realization. No model refitting is required; the difference between LR and CB lies only in how valid parameter realizations are generated before the common propagation procedure is applied.

\paragraph{Bootstrap methods}

Bootstrap methods quantify epistemic uncertainty by repeatedly generating bootstrap datasets, refitting the corresponding statistical model, and constructing confidence intervals and confidence bands from the resulting bootstrap distributions. Unlike analytical approximations and Monte Carlo sampling from the asymptotic parameter distribution, bootstrap methods repeatedly re-estimate the model parameters and therefore naturally account for the variability introduced by the model-fitting procedure itself. Two bootstrap methods were considered: nonparametric bootstrap (NPBS) and parametric bootstrap (PBS) \citep{Efron1994}.

Let

\begin{equation}
\mathcal{D}
=
\left\{
(x_i,y_i)
\right\}_{i=1}^{n},
\label{eq:bootstrap_dataset}
\end{equation}

denote the observed dataset, where $x_i$ is the biomarker value, $y_i\in\{0,1\}$ is the corresponding binary outcome for patient $i$, and $n$ is the total number of patients. Let $b=1,\ldots,B$ denote the bootstrap replicate index, where $B$ is the total number of bootstrap replicates. The bootstrap dataset corresponding to replicate $b$ is denoted by

\begin{equation}
\mathcal{D}^{(b)},
\label{eq:bootstrap_replicate}
\end{equation}

whose construction depends on the selected bootstrap method and is described below.

For each bootstrap replicate, the corresponding statistical model, $m\in\{\mathrm{LR},\mathrm{CB}\}$, was refitted to obtain the bootstrap parameter estimate $\hat{\boldsymbol{\Theta}}_{m}^{(b)}$, from which the estimated risk function $p(x;\hat{\boldsymbol{\Theta}}_{m}^{(b)})$ and the corresponding model characteristics $x_{50,m}^{(b)}$ and $s_{50,m}^{(b)}$ were derived.

Confidence intervals and confidence bands were constructed using the percentile bootstrap method. For an arbitrary estimated quantity $\psi$, let

\[
\psi^{(1)},
\psi^{(2)},
\ldots,
\psi^{(B)}
\]

denote its bootstrap estimates obtained from the $B$ bootstrap replicates. The corresponding two-sided $(1-\alpha)\times100\%$ confidence interval was defined as

\begin{equation}
CI_{1-\alpha}(\psi)
=
\left[
Q_{\alpha/2}(\psi),
Q_{1-\alpha/2}(\psi)
\right],
\label{eq:bootstrap_percentile}
\end{equation}

where $Q_p(\psi)$ denotes the empirical $p$-quantile of the bootstrap distribution of $\psi$. Throughout this study, $\alpha=0.05$, corresponding to 95\% confidence intervals.

\paragraph{Nonparametric bootstrap (NPBS)}

In the nonparametric bootstrap (NPBS), bootstrap datasets were generated by sampling complete observations from the observed dataset with replacement, without assuming any parametric form for the underlying data-generating process \citep{EfronTibshirani1994}. Thus, for each bootstrap replicate $b=1,\ldots,B$, the bootstrap dataset

\begin{equation}
\mathcal{D}^{(b)}
=
\left\{
(\tilde{x}_i^{(b)},\tilde{y}_i^{(b)})
\right\}_{i=1}^{n},
\label{eq:npbs_dataset}
\end{equation}

consisted of $n$ observation pairs sampled with replacement from the original dataset $\mathcal{D}$, where $(\tilde{x}_i^{(b)},\tilde{y}_i^{(b)})$ denotes the $i$-th resampled observation in the $b$-th bootstrap dataset. Bootstrap datasets containing only one outcome class were discarded because the corresponding statistical model could not be estimated. Additional bootstrap datasets were generated until the required number of bootstrap replicates was obtained.

For each retained bootstrap dataset $\mathcal{D}^{(b)}$, both model formulations were refitted using their respective estimation procedures. The LR model was refitted by maximum-likelihood estimation, whereas the CB model was refitted by maximizing the posterior objective function over the feasible parameter region. The resulting bootstrap estimates $\hat{\boldsymbol{\Theta}}_{\mathrm{LR}}^{(b)}$ and $\hat{\boldsymbol{\Theta}}_{\mathrm{CB}}^{(b)}$ were used to calculate the corresponding fitted risk functions and derived model characteristics.

\paragraph{Parametric bootstrap (PBS)}

In the parametric bootstrap, bootstrap datasets were generated from the fitted probabilistic model rather than by resampling observations from the original dataset \citep{EfronTibshirani1994}. For each bootstrap replicate $b=1,\ldots,B$, a synthetic dataset was generated using the parameter estimates obtained from the original data, and the corresponding model was subsequently refitted.

The construction of the synthetic dataset reflects the probabilistic formulation of each model. For LR, the observed biomarker values were kept fixed and binary outcomes were generated from the fitted risk function according to

\begin{equation}
Y_{i,\mathrm{LR}}^{(b)}
\sim
\mathrm{Bernoulli}
\!\left(
p(x_i;\hat{\boldsymbol{\Theta}}_{\mathrm{LR}})
\right),
\qquad
i=1,\ldots,n,
\label{eq:pbs_lr_generation}
\end{equation}

giving the bootstrap dataset

\begin{equation}
\mathcal{D}_{\mathrm{LR}}^{(b)}
=
\left\{
(x_i,Y_{i,\mathrm{LR}}^{(b)})
\right\}_{i=1}^{n}.
\label{eq:pbs_lr_dataset}
\end{equation}

For the CB model, both the binary outcome and biomarker value were generated from the fitted generative model. The outcome was first sampled according to

\begin{equation}
Y_{i,\mathrm{CB}}^{(b)}
\sim
\mathrm{Bernoulli}
\left(
\hat{\omega}
\right),
\qquad
i=1,\ldots,n,
\label{eq:pbs_cb_label}
\end{equation}

where $\hat{\omega}$ is the adverse-event prevalence estimated from the fitted CB model. Conditional on the generated outcome, the corresponding biomarker value was sampled from the fitted group-specific distribution,

\begin{equation}
X_{i,\mathrm{CB}}^{(b)}
\mid
Y_{i,\mathrm{CB}}^{(b)}
\sim
\begin{cases}
\Gamma(\hat{k},\hat{\vartheta}),
&
Y_{i,\mathrm{CB}}^{(b)}=0,
\\[2mm]
\mathrm{BetaPrime}(\hat{a},\hat{b},\hat{s}),
&
Y_{i,\mathrm{CB}}^{(b)}=1,
\end{cases}
\label{eq:pbs_cb_generation}
\end{equation}

giving the bootstrap dataset

\begin{equation}
\mathcal{D}_{\mathrm{CB}}^{(b)}
=
\left\{
(X_{i,\mathrm{CB}}^{(b)},Y_{i,\mathrm{CB}}^{(b)})
\right\}_{i=1}^{n}.
\label{eq:pbs_cb_dataset}
\end{equation}

For both model formulations, bootstrap datasets containing only one outcome class were discarded and additional datasets were generated until the required number of bootstrap replicates was obtained. Each retained dataset was then refitted using the same estimation procedure applied to the original data: maximum-likelihood estimation for LR and posterior maximization over the feasible parameter region for CB. The resulting bootstrap parameter estimates were propagated to the corresponding fitted risk functions and derived characteristics. The empirical bootstrap distributions were then used to construct the percentile confidence intervals and pointwise confidence bands defined above.

% -------------------------------------------------------------------------
\subsubsection{Noise-induced uncertainty}

Measurement-noise propagation was evaluated to examine how measurement variability in the $^{18}$F-FDG PET/CT-derived quantitative imaging biomarker propagates to the estimated risk. The perturbation magnitudes were selected based on a preceding characterization of noise in the PET data. In that analysis, an approximately homogeneous region defined by $1.0<\mathrm{SUV}<1.2$ was used to characterize image noise. The additive noise level was characterized by the standard deviation of the SUV values within this region, whereas the multiplicative noise level was characterized by the corresponding coefficient of variation. This analysis yielded representative perturbation magnitudes of $\sigma_{\mathrm{add}}=0.144$ and $\sigma_{\mathrm{mult}}=0.129$, which were subsequently used as fixed inputs in the uncertainty-propagation analysis.

Noise was propagated at the level of the subject-level biomarker $X$, defined as the maximum $\mathrm{SUV}_{95\%}$ across the available longitudinal visits. Measurement variability was represented using additive and multiplicative Gaussian perturbations. For the additive model, the perturbed biomarker value for subject $i$ in realization $m$ was generated as
%
\begin{equation}
\tilde{x}_i^{(m)}
=
x_i+\varepsilon_i^{(m)},
\qquad
\varepsilon_i^{(m)}
\sim
\mathcal{N}(0,\sigma_{\mathrm{add}}^2),
\label{eq:meas_add}
\end{equation}
%
whereas for the multiplicative model,
%
\begin{equation}
\tilde{x}_i^{(m)}
=
x_i\exp\!\left(\varepsilon_i^{(m)}\right),
\qquad
\varepsilon_i^{(m)}
\sim
\mathcal{N}(0,\sigma_{\mathrm{mult}}^2).
\label{eq:meas_mult}
\end{equation}
%
Positive perturbed biomarker values were retained by applying a small positive lower bound when required. Two complementary analyses were performed to examine how these perturbations propagated to the predicted risk: a fixed-model analysis and a model-refitting analysis.

\paragraph{Fixed-model analysis.}

The fixed-model analysis quantified the effect of biomarker measurement variability on predicted risk while keeping the model parameters estimated from the unperturbed data unchanged. For each Monte Carlo realization, biomarker values along the evaluation grid were perturbed according to the selected additive or multiplicative noise model and passed through the fitted risk function. Thus, variation in the resulting predicted risk reflects propagation of biomarker measurement variability through a fixed fitted model.

At each biomarker value, the resulting predicted probabilities across perturbation realizations were summarized using the 2.5th, 50th, and 97.5th percentiles, with the 2.5th and 97.5th percentiles defining the 95\% uncertainty band.

\paragraph{Model-refitting analysis.}

The model-refitting analysis evaluated the effect of biomarker measurement variability when model estimation was repeated after perturbation. For each Monte Carlo realization, the observed subject-level biomarker values were independently perturbed according to the selected noise model, while the corresponding outcome labels were kept unchanged. The statistical model was then refitted to the perturbed dataset,
%
\[
\left\{
\left(\tilde{x}_i^{(m)},y_i\right)
\right\}_{i=1}^{n},
\]
%
and the resulting fitted risk function was evaluated over the original biomarker grid. At each biomarker value, the 2.5th, 50th, and 97.5th percentiles of the predicted probabilities across successful refits were obtained, with the 2.5th and 97.5th percentiles defining the 95\% uncertainty band.

For comparison of noise propagation across models, datasets, and perturbation types, the vertical width of each 95\% uncertainty band was evaluated at the midpoint $x_{50}$ of the risk function fitted to the unperturbed data, defined by
%
\[
p(x_{50};\hat{\bTheta})=0.5.
\]
%
The midpoint was held fixed for this comparison. The bandwidth, denoted by $\Delta$, was calculated as the difference between the 97.5th and 2.5th percentiles of the propagated predicted probabilities evaluated at $x_{50}$. Thus, $\Delta$ quantifies variability in predicted risk at the clean-data midpoint rather than uncertainty in $x_{50}$ itself.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}

The probabilistic risk-modeling and uncertainty-quantification framework was evaluated using the FULL and TRIM datasets of patients receiving immune checkpoint inhibitor therapy. The fitted risk models are presented first, followed by model evaluation, derived risk-function characteristics, sensitivity analysis, and uncertainty quantification.

\subsection{Model fitting}

\subsubsection{Standard and cubic logistic regression models}

The fitted risk functions obtained using the standard and cubic logistic regression models are shown in Figure~\ref{fig:logistic_fits}.

For both modeling approaches, the predicted probability of pneumonitis increased with increasing biomarker values. Differences between the FULL and TRIM datasets were observed on both the log-transformed and raw biomarker scales. Removal of the identified normal-control observation resulted in changes in the fitted risk functions for both modeling approaches.

The standard and cubic logistic regression models produced similar overall risk patterns, although differences were observed over the biomarker range where the predicted probability changed most rapidly. These differences were more apparent in the FULL dataset and were reduced in the TRIM dataset. At lower and higher biomarker values, the fitted models showed more similar behavior.

\begin{figure}[!htbp]
\centering

\includegraphics[width=\figWmed]{Figures/LOGISTICFITS.pdf}

\vspace{0.4cm}

\includegraphics[width=\figWmed]{Figures/cubic_logistic_fits.pdf}

\caption{
Standard logistic regression (Panels A and B) and cubic logistic regression (Panels C and D) fits for the FULL and TRIM datasets.
The conditional AE probability,
\(p(x)=\Prob(\AEgrp\mid X=x)\),
is shown as a function of the lung biomarker \(x\), defined as the maximum longitudinal \(\SUV_{95\%}\) value across all eligible visits.
Panels A and C show the log-transformed biomarker scale, whereas Panels B and D show the raw biomarker scale.
Blue and orange curves correspond to the FULL and TRIM datasets, respectively. The TRIM dataset excludes the identified normal-control observation.
Observed binary outcomes are shown as points. A small vertical jitter was added for visualization only.
}
\label{fig:logistic_fits}
\end{figure}

To examine the local behavior of the cubic logistic objective function around the fitted solution, profiles of the negative log-likelihood are shown in Figure~\ref{fig:cubic_logistic_cost_min}. The profiles show differences in local curvature across the model parameters. Some parameter directions exhibit relatively sharp minima, whereas others are flatter or less regular. These differences indicate that the local behavior of the objective function varies across parameter directions and are consistent with differences in the precision with which individual parameters are estimated.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=\figWmed]{logit_cost_min_paper.pdf}
  \caption{Negative log-likelihood profiles around the fitted solution for the cubic logistic model. Each curve shows the change in the objective function when a single parameter in the constrained parameterization is perturbed while the remaining parameters are held fixed. Panel (A) shows the FULL dataset and panel (B) the TRIM dataset using the log-transformed biomarker. The vertical axis is shown on a logarithmic scale.}
  \label{fig:cubic_logistic_cost_min}
\end{figure}
\FloatBarrier
% -------------------------------------------------------------------------
\subsubsection{Constrained Bayesian model}

The fitted risk functions obtained using the constrained Bayesian model are shown in Figure~\ref{fig:bayes_fits}.

For both datasets, the predicted probability of pneumonitis increased with increasing biomarker values. Differences were observed between the risk functions fitted to the FULL and TRIM datasets, particularly over the biomarker range where the predicted probability changed most rapidly. At lower and higher biomarker values, the fitted risk functions showed more similar behavior.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.5\figWfull]{BAYESIANFITs.pdf}
\caption{
Constrained Bayesian model fits for the FULL and TRIM datasets.
The conditional AE probability,
\(p(x)=\Prob(\AEgrp\mid X=x)\),
is shown as a function of the lung biomarker \(x\), defined as the maximum longitudinal \(\SUV_{95\%}\) value across all eligible visits.
Blue and orange curves correspond to the FULL and TRIM datasets, respectively. The TRIM dataset excludes the identified normal-control observation.
Observed binary outcomes are overlaid as points. Filled markers denote the FULL dataset, and outlined markers denote the TRIM dataset. A small vertical jitter was added for visualization only.
}
  \label{fig:bayes_fits}
\end{figure}
\FloatBarrier

% -------------------------------------------------------------------------
\subsection{Model evaluation}

The model evaluation measures are summarized in Table~\ref{tab:gof}.

\begin{table}[!htbp]
\centering
\caption{
Model evaluation using the Bernoulli log-likelihood (LLF), Akaike information criterion (AIC), and Bayesian information criterion (BIC). Larger LLF values indicate better agreement with the observed outcomes, whereas smaller AIC and BIC values indicate a more favorable balance between fit and model complexity.
}
\label{tab:gof}

\begin{tabular}{llcccc}
\toprule
Model & Dataset & $q$ & LLF $\uparrow$ & AIC $\downarrow$ & BIC $\downarrow$\\
\midrule
Logistic (RAW) & FULL & 2 & $-7.394$ & 18.788 & 22.909\\
Logistic (LOG) & FULL & 2 & $-6.795$ & 17.591 & 21.712\\
Logistic (RAW) & TRIM & 2 & $-2.769$ & 9.539 & 13.625\\
Logistic (LOG) & TRIM & 2 & $-2.712$ & 9.424 & 13.510\\
Constrained Bayesian & FULL & 6 & $-4.687$ & 21.373 & 33.736\\
Constrained Bayesian & TRIM & 6 & $-2.610$ & 17.220 & 29.478\\
\bottomrule
\end{tabular}
\end{table}

Among the standard logistic regression models, the log-transformed formulation yielded slightly higher LLF values and lower AIC and BIC values than the corresponding raw-scale formulation for both the FULL and TRIM datasets. The differences between the two logistic formulations were smaller in the TRIM dataset. Compared with the FULL dataset, the TRIM dataset yielded higher LLF values and lower AIC and BIC values for both logistic formulations.

The constrained Bayesian model yielded higher LLF values than the standard logistic regression models for both datasets. However, its AIC and BIC values were also higher because of the larger number of estimated parameters.

Among the models for which AIC and BIC were reported, the log-transformed logistic regression model had the lowest AIC and BIC values for both datasets, whereas the constrained Bayesian model had the highest LLF values.

% -------------------------------------------------------------------------
\subsection{Risk function characteristics}

The fitted risk functions were characterized using the midpoint \(x_{50}\) and midpoint slope \(s_{50}\), defined in Eqs.~\eqref{eq:x50} and \eqref{eq:s50}, respectively. These quantities provide common summaries of the fitted risk functions that can be compared across model formulations without direct comparison of model-specific parameters. The resulting values are reported in Table~\ref{tab:curve_characteristics}.

The midpoint \(x_{50}\) is the biomarker value corresponding to a predicted risk of 50\%, whereas \(s_{50}\) is the rate of change in predicted risk with respect to the biomarker at \(x_{50}\). The midpoint is used here as a characteristic of the fitted risk function and should not be interpreted as a validated clinical decision threshold.

\begin{table}[!htbp]
\centering
\caption{
Risk-function characteristics for the standard logistic regression and constrained Bayesian models. The midpoint \(x_{50}\) and midpoint slope \(s_{50}\) were computed according to Eqs.~\eqref{eq:x50} and \eqref{eq:s50}, respectively.
}
\label{tab:curve_characteristics}
\begin{tabular}{ll
                S[table-format=1.3]
                S[table-format=1.3]}
\toprule
Dataset & Model & $x_{50}$ & $s_{50}$ \\
\midrule
FULL & Logistic (raw \(x\)) & 2.054 & 1.332 \\
FULL & Logistic (log \(x\)) & 2.016 & 2.684 \\
\midrule
TRIM & Logistic (raw \(x\)) & 1.775 & 3.512 \\
TRIM & Logistic (log \(x\)) & 1.769 & 6.213 \\
\midrule
FULL & Constrained Bayesian & 1.792 & 1.830 \\
TRIM & Constrained Bayesian & 1.754 & 3.356 \\
\bottomrule
\end{tabular}
\end{table}

Across the standard logistic regression models, the estimated midpoint \(x_{50}\) was 2.054 and 2.016 for the raw-scale and log-transformed formulations, respectively, in the FULL dataset, and 1.775 and 1.769 in the TRIM dataset (Table~\ref{tab:curve_characteristics}). The two logistic formulations therefore produced similar \(x_{50}\) values within each dataset. In contrast, their \(s_{50}\) values differed, with larger values for the log-transformed formulation in both datasets.

The constrained Bayesian model yielded \(x_{50}\) values of 1.792 and 1.754 for the FULL and TRIM datasets, respectively, while the corresponding \(s_{50}\) values were 1.830 and 3.356. Across the reported model formulations and datasets, the numerical variation in \(x_{50}\) was smaller than that observed for \(s_{50}\).
% -------------------------------------------------------------------------
\subsubsection{Sensitivity analysis of risk-function characteristics}
\label{sec:log_sens_results}

Sensitivity and elasticity analyses were performed as described in Section~\ref{sec:sensitivity_elasticity}. The results are presented separately for the standard logistic regression and constrained Bayesian models.

\paragraph{Logistic regression.}

The absolute elasticities,
\(|\mathcal{E}_{\theta_j}(x_{50})|\) and
\(|\mathcal{E}_{\theta_j}(s_{50})|\),
with respect to the intercept and slope parameters
\((\theta_0,\theta_1)\)
for the four logistic model variants
(FULL/TRIM \(\times\) RAW/LOG)
are shown in Figure~\ref{fig:elas_logistic}.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\figWwide]{SENLOGISTIC.pdf}
\caption{
Elasticity analysis of the standard logistic regression models.
Panel~(A) shows the absolute elasticities
\(|\mathcal{E}_{\theta_j}(x_{50})|\),
and Panel~(B) shows
\(|\mathcal{E}_{\theta_j}(s_{50})|\),
with respect to the parameter vector
\(\bTheta=(\theta_0,\theta_1)\),
where \(\theta_0\) denotes the intercept and
\(\theta_1\) the slope parameter.
Bars represent the FULL and TRIM datasets under the raw-scale and log-transformed formulations.
}
\label{fig:elas_logistic}
\end{figure}

For the midpoint \(x_{50}\), the absolute elasticities with respect to \(\theta_0\) and \(\theta_1\) were equal within each model variant. The corresponding values were approximately \(1\) for the raw-scale models and approximately \(0.57\)--\(0.70\) for the log-transformed models. Differences between the FULL and TRIM datasets were comparatively small.

For the midpoint slope \(s_{50}\), a different pattern was observed. In the raw-scale models, the absolute elasticity with respect to \(\theta_0\) was zero, whereas that with respect to \(\theta_1\) was approximately \(2\). In the log-transformed models, both parameters had non-zero absolute elasticities, with values of approximately \(0.57\)--\(0.70\) for \(\theta_0\) and \(1.57\)--\(1.70\) for \(\theta_1\). Differences between the FULL and TRIM datasets were smaller than those between the raw-scale and log-transformed formulations.

\paragraph{Constrained Bayesian model.}

The absolute elasticities,
\(|\mathcal{E}_{\theta_j}(x_{50})|\) and
\(|\mathcal{E}_{\theta_j}(s_{50})|\),
with respect to the parameter vector
\(\bTheta=(\omega,a,b,s,k,\vartheta)\)
for the FULL and TRIM datasets are shown in Figure~\ref{fig:elas_bayes}.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\figWwide]{SENBAYES.pdf}
\caption{
Elasticity analysis of the constrained Bayesian model.
Panel~(A) shows the absolute elasticities
\(|\mathcal{E}_{\theta_j}(x_{50})|\),
and Panel~(B) shows
\(|\mathcal{E}_{\theta_j}(s_{50})|\),
with respect to the parameter vector
\(\bTheta=(\omega,a,b,s,k,\vartheta)\).
The vertical axis is displayed on a logarithmic scale.
Bars represent the FULL and TRIM datasets.
}
\label{fig:elas_bayes}
\end{figure}

For both \(x_{50}\) and \(s_{50}\), the smallest absolute elasticities were associated with the prevalence parameter \(\omega\), whereas the largest were observed for the Gamma-distribution parameters \(k\) and \(\vartheta\). The parameters \(a\), \(b\), and \(s\) showed intermediate absolute elasticities.

For the parameters shown, the FULL dataset yielded larger absolute elasticities than the TRIM dataset. Absolute elasticities were also generally larger for \(s_{50}\) than for \(x_{50}\), indicating a greater relative response of \(s_{50}\) to local proportional changes in the fitted model parameters.
% -------------------------------------------------------------------------
\subsection{Uncertainty Quantification}
\subsubsection{Confidence intervals and confidence bands}

Uncertainty in the logistic regression model was evaluated using four complementary approaches: the Wald approximation, Monte Carlo approximation (MCA), the nonparametric pairs bootstrap (NPBS), and the parametric bootstrap (PBS). The Wald approximation estimates uncertainty analytically. MCA estimates uncertainty by propagating the estimated parameter distribution through the fitted risk function and its derived risk function characteristics. NPBS estimates uncertainty empirically by repeated resampling of the observed data, whereas PBS estimates uncertainty empirically by repeated sampling from the fitted probabilistic model.

Uncertainty is summarized using confidence bands for the fitted risk function and confidence intervals for its derived risk function characteristics. These summaries describe different but complementary aspects of model uncertainty. Confidence bands characterize uncertainty in the fitted risk function across the biomarker range, whereas confidence intervals for \(x_{50}\) characterize uncertainty in the estimated transition location of the fitted risk function.

Because \(x_{50}\) is the primary clinically interpretable characteristic of the fitted risk function in this study, confidence intervals for \(x_{50}\) provide the principal inferential summary. Confidence bands are used primarily to visualize uncertainty in the fitted risk function and to assess how uncertainty varies across the biomarker range.

\paragraph{Logistic regression.}

The bootstrap diagnostics are shown in Figure~\ref{fig:bootstraps}. The NPBS parameter distribution exhibits multiple modes (Figure~\ref{fig:bootstraps}a). Focusing on the dominant mode, the maximum-likelihood estimate lies near the highest-density region, and the local Gaussian approximation around \(\hat{\boldsymbol{\Theta}}\) captures the location and spread of this dominant mode (Figure~\ref{fig:bootstraps}b).

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[t]{0.49\textwidth}
    \centering
    \includegraphics[width=\linewidth]{logit_nonpar_boots.pdf}
    \label{fig:bootstraps:a}
  \end{minipage}
  \hfill
  \begin{minipage}[t]{0.49\textwidth}
    \centering
    \includegraphics[width=\linewidth]{logit_nonpar_boots_sel.pdf}
    \label{fig:bootstraps:b}
  \end{minipage}
  \caption{
Heat maps of the NPBS distribution of the logistic regression parameter estimates \((\theta_0,\theta_1)\), obtained by maximum likelihood estimation (MLE) for the FULL cohort using the RAW predictor formulation. Panel~(a) shows the distribution over the complete parameter range. In contrast, Panel~(b) provides a magnified view of the dominant mode around the MLE obtained from the original data, indicated by the red point. In both panels, the upper and right subpanels show the marginal distributions of \(\theta_0\) and \(\theta_1\), respectively. In Panel~( b), the red curves show the marginal densities implied by the local Gaussian approximation of the parameter estimates; see Eq.~\eqref{eq:mc_sampling}.
}
  \label{fig:bootstraps}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=\figWfull]{LOGISTICCIs.pdf}
\caption{
Fitted logistic regression risk functions with 95\% confidence bands for the FULL and TRIM datasets. Panels~(A)--(B) show the FULL dataset under the LOG and RAW formulations, respectively, and Panels~(C)--(D) show the corresponding TRIM analyses. The vertical axis represents the predicted probability of an adverse event, \(p(x)=P(\mathrm{AE}\mid X=x)\). The biomarker \(X\) is defined for each subject as the maximum across all longitudinal visits of \(\mathrm{SUV}_{95\%}\), the 95th percentile of \(^{18}\mathrm{F}\)-FDG uptake within the lung region; the horizontal axis shows \(\log(X)\) for the LOG formulation and \(X\) for the RAW formulation. Solid curves represent the fitted risk functions. Confidence bands were obtained using the Wald approximation, Monte Carlo approximation (MCA), nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS), as identified by the colours and line styles. Observed normal-control (NC) and adverse-event (AE) data are overlaid with a small vertical jitter to reduce overlap.
}
  \label{fig:logistic_cis}
\end{figure}
Coefficient estimates and corresponding 95\% confidence intervals are summarized in Table~\ref{tab:logistic_param_cis}. In the FULL dataset, the NPBS and PBS intervals are substantially wider than the Wald and MCA intervals for both regression coefficients. In the TRIM dataset, the NPBS and PBS intervals are narrower than their corresponding intervals in the FULL dataset.

\begin{table}[!htbp]
\centering
\caption{Estimated logistic regression parameters with corresponding 95\% confidence intervals obtained using the Wald approximation, Monte Carlo approximation (MCA), nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS).}
\label{tab:logistic_param_cis}

\resizebox{\textwidth}{!}{%
\begin{tabular}{llccccccc}
\toprule
Panel & Parameter & Estimate & Wald & MCA & NPBS & PBS \\
\midrule

A & $\theta_0$ & -7.527 & [-11.932, -3.123] & [-11.924, -3.137] & [-134.869, -5.297] & [-1350.370, -4.714] \\
  & $\theta_1$ & 10.738 & [3.403, 18.073] & [3.425, 18.067] & [5.849, 281.074] & [5.481, 1821.060] \\

\addlinespace

B & $\theta_0$ & -11.684 & [-18.564, -4.803] & [-18.550, -4.837] & [-279.442, -7.344] & [-1888.580, -7.116] \\
  & $\theta_1$ & 5.687 & [1.857, 9.517] & [1.867, 9.508] & [2.855, 172.783] & [2.985, 899.685] \\

\addlinespace

C & $\theta_0$ & -5.741 & [-8.270, -3.211] & [-8.270, -3.224] & [-6.341, -5.070] & [-5.927, -3.348] \\
  & $\theta_1$ & 8.288 & [3.459, 13.118] & [3.475, 13.109] & [4.995, 9.921] & [3.277, 8.288] \\

\addlinespace

D & $\theta_0$ & -9.477 & [-13.988, -4.966] & [-13.975, -4.991] & [-10.699, -7.223] & [-9.591, -4.646] \\
  & $\theta_1$ & 4.715 & [1.942, 7.488] & [1.947, 7.484] & [2.412, 5.683] & [1.603, 4.778] \\

\bottomrule
\end{tabular}%
}
\end{table}

Midpoint estimates \(x_{50}\) and their corresponding 95\% confidence intervals are summarized in Table~\ref{tab:logistic_x50_cis_all}. These intervals quantify uncertainty in the biomarker location at which the fitted risk reaches 50\%, rather than vertical uncertainty in the predicted probability at \(x=x_{50}\). The point estimates are similar for the RAW and LOG formulations within each dataset, with values of approximately 2.0 SUV for the FULL dataset and 1.77 SUV for the TRIM dataset. The confidence-interval widths vary across uncertainty quantification methods and between the FULL and TRIM datasets.


\begin{table}[!htbp]
\centering
\caption{Midpoint estimates $\widehat{x}_{50}$ (SUV units) with corresponding 95\% confidence intervals obtained using the Wald approximation, Monte Carlo approximation (MCA), the nonparametric pairs bootstrap (NPBS), and the parametric bootstrap (PBS).}
\label{tab:logistic_x50_cis_all}

\begin{tabular}{llccccc}
\toprule
Dataset & Predictor & $\widehat{x}_{50}$ & Wald & MCA & NPBS & PBS \\
\midrule
FULL & RAW & 2.055 & [1.735, 2.374] & [1.804, 2.767] & [1.612, 2.798] & [1.730, 2.495] \\
FULL & LOG & 2.016 & [1.711, 2.375] & [1.759, 2.872] & [1.611, 2.775] & [1.709, 2.508] \\
\midrule
TRIM & RAW & 1.775 & [1.571, 1.979] & [1.613, 2.576] & [1.607, 2.002] & [1.618, 1.936] \\
TRIM & LOG & 1.769 & [1.579, 1.981] & [1.602, 2.651] & [1.605, 1.976] & [1.617, 1.930] \\
\bottomrule
\end{tabular}
\end{table}


\paragraph{Monotone cubic logistic extension.}

Figure~\ref{fig:logit_cubic_ci} shows the 95\% confidence bands for the monotone cubic logistic model obtained using the Wald approximation, Monte Carlo approximation (MCA), the nonparametric pairs bootstrap (NPBS), and the parametric bootstrap (PBS). The confidence-band widths and shapes differed among the four uncertainty quantification methods. These differences were more pronounced in the FULL dataset, whereas the confidence bands obtained using the different methods were more similar in the TRIM dataset.


\begin{figure}[!htbp]
  \centering
  \includegraphics[width=\figWfull]{Figures/logit_fit_CI_paper.pdf}
\caption{Fitted cubic logistic risk functions with 95\% confidence bands for the FULL and TRIM datasets. Panels~(A) and~(C) show the FULL and TRIM datasets, respectively, under the LOG predictor formulation, whereas Panels~(B) and~(D) show the corresponding analyses under the RAW predictor formulation. The vertical axis represents the predicted probability of an adverse event, \(p(x)=P(\mathrm{AE}\mid X=x)\). The biomarker \(X\) is defined for each subject as the maximum across all longitudinal visits of \(\mathrm{SUV}_{95\%}\), the 95th percentile of \(^{18}\mathrm{F}\)-FDG uptake within the lung region; the horizontal axis shows \(\log(X)\) for the LOG formulation and \(X\) for the RAW formulation. Solid curves represent the fitted risk functions. Confidence bands were obtained using the Wald approximation, Monte Carlo approximation (MCA), nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS), as identified by the colors and line styles. Observed normal-control (NC) and adverse-event (AE) data are overlaid with a small vertical jitter to reduce overlap.}
  \label{fig:logit_cubic_ci}
\end{figure}

\paragraph{Constrained Bayesian model.}

For the constrained Bayesian model, uncertainty was quantified around the maximum a posteriori (MAP) estimate of the fitted risk function.

\[
\hat p(x)=p(x;\hat{\bTheta}),
\]

using three uncertainty quantification approaches: a Wald approximation, a nonparametric pairs bootstrap (NPBS), and a parametric bootstrap (PBS). As in the logistic regression analysis, uncertainty was summarized at three levels: interval estimates for the model parameters, confidence bands for the fitted risk function \(p(x)\), and interval estimates for the characteristic point \(x_{50}\).

Parameter interval estimates obtained using the three uncertainty quantification approaches are presented in Table~\ref{tab:cb_param_ci}. The interval estimates vary across model parameters and among the three uncertainty quantification approaches. Wide interval estimates are observed for several parameters in both the FULL and TRIM datasets.

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{5pt}
\renewcommand{\arraystretch}{1.08}
\caption{Parameter estimates and corresponding 95\% confidence intervals obtained using the Wald approximation, nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS) for the FULL and TRIM datasets.}
\label{tab:cb_param_ci}
\begin{tabular}{@{}c c S[table-format=3.4] l l l@{}}
\toprule
{Dataset} & {Parameter} & {{Estimate}} & {Wald (95\% CI)} & {NPBS (95\% CI)} & {PBS (95\% CI)} \\
\midrule
FULL & $a$ & 183.32 & [21.58, 549.25] & [32.16, 493.64] & [25.47, 436.25] \\
FULL & $b$ & 12.74 & [0.0630, 27.99] & [7.10, 518.05] & [6.99, 443.62] \\
FULL & $s$ & 0.1560 & [0.0117, 1.2015] & [0.0853, 16.05] & [0.0701, 20.05] \\
FULL & $k$ & 27.81 & [17.13, 38.08] & [16.57, 56.13] & [20.03, 42.93] \\
FULL & $\vartheta$ & 0.0470 & [0.0063, 0.0634] & [0.0235, 0.0805] & [0.0303, 0.0653] \\
\addlinespace[2pt]
TRIM & $a$ & 207.61 & [26.34, 1745.28] & [41.12, 436.72] & [34.90, 391.87] \\
TRIM & $b$ & 13.06 & [0.0452, 29.85] & [7.86, 626.35] & [8.07, 551.22] \\
TRIM & $s$ & 0.1410 & [3.42e-05, 6.80] & [0.0755, 17.95] & [0.0688, 10.63] \\
TRIM & $k$ & 37.55 & [23.18, 51.99] & [26.91, 59.86] & [27.17, 58.89] \\
TRIM & $\vartheta$ & 0.0342 & [4.11e-05, 0.1113] & [0.0217, 0.0467] & [0.0220, 0.0475] \\
\bottomrule
\end{tabular}
\end{table}

Interval estimates for the characteristic point \(x_{50}\) are presented in Table~\ref{tab:cb_x50_ci}. The estimated \(x_{50}\) values are close across the three uncertainty quantification approaches for both datasets, whereas the corresponding interval widths differ.

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{5pt}
\renewcommand{\arraystretch}{1.08}
\caption{Estimated \(x_{50}\) values and corresponding 95\% confidence intervals obtained using the Wald approximation, nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS) for the FULL and TRIM datasets.}
\label{tab:cb_x50_ci}
\begin{tabular}{@{}c S[table-format=1.3] l l l@{}}
\toprule
{Dataset} & {{$\hat x_{50}$}} & {Wald (95\% CI)} & {NPBS (95\% CI)} & {PBS (95\% CI)} \\
\midrule
FULL & 1.792 & [1.647, 2.919] & [1.609, 2.888] & [1.669, 2.643] \\
TRIM & 1.754 & [1.638, 1.936] & [1.608, 2.009] & [1.651, 3.478] \\
\bottomrule
\end{tabular}
\end{table}

Figure~\ref{fig:BayesianCIs} shows the corresponding 95\% confidence bands for the fitted risk function. For all three uncertainty quantification approaches, the confidence bands are widest over the range where the fitted probability changes most rapidly. Outside this range, the confidence bands narrow. The confidence bands for the TRIM dataset are generally narrower than those for the FULL dataset.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\figWfull]{BayesianCIs.pdf}
\caption{
Fitted constrained Bayesian risk functions with corresponding 95\% confidence bands for the FULL (A) and TRIM (B) datasets. The vertical axis represents the predicted probability of an adverse event, \(p(x)=P(\mathrm{AE}\mid X=x)\), and the horizontal axis represents the biomarker \(X\), defined for each subject as the maximum across all longitudinal visits of \(\mathrm{SUV}_{95\%}\), where \(\mathrm{SUV}_{95\%}\) is the 95th percentile of \(^{18}\mathrm{F}\)-FDG uptake within the lung region. The solid curve represents the fitted risk function. Shaded regions show the 95\% confidence bands obtained using the Wald approximation, nonparametric pairs bootstrap (NPBS), and parametric bootstrap (PBS). Observed normal-control (NC) and adverse-event (AE) data are overlaid with a small vertical jitter to reduce overlap.
}
\label{fig:BayesianCIs}
\end{figure}

%-----------------------------------------------------------------------------------------------------------------
\subsubsection{Noise-induced uncertainty}
\label{sec:uq_results:noise}

Measurement-noise propagation was evaluated for the logistic regression and constrained Bayesian models using the additive and multiplicative perturbations described in the Methods, with $\sigma_{\mathrm{add}}=0.144$ and $\sigma_{\mathrm{mult}}=0.129$, respectively. For each model, dataset, and noise type, 95\% uncertainty bands were obtained from the fixed-model and model-refitting analyses. The uncertainty-band width, denoted by $\Delta$, was evaluated at the clean-data midpoint $x_{50}$. This quantity represents vertical variability in predicted probability at the fixed midpoint and should therefore be distinguished from uncertainty in the estimated location of $x_{50}$.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\figWfull]{Figures/NOISELOGISTIC.pdf}
\caption{Measurement-noise propagation through the logistic regression model for the FULL (A--B) and TRIM (C--D) datasets. Panels A and C show multiplicative Gaussian perturbations ($\sigma_{\mathrm{mult}}=0.129$), and Panels B and D show additive Gaussian perturbations ($\sigma_{\mathrm{add}}=0.144$). The black curve represents the model fitted to the unperturbed data. Dark and light shaded regions represent the 95\% uncertainty bands from the model-refitting and fixed-model analyses, respectively. The vertical dashed line indicates the clean-data midpoint $x_{50}$. Observed NC and AE data are shown at probabilities 0 and 1, respectively.}
\label{fig:noise_logistic}
\end{figure}

For the logistic regression model (Figure~\ref{fig:noise_logistic}), the clean-data midpoint was $x_{50}=2.02$ for the FULL dataset and $x_{50}=1.77$ for the TRIM dataset. In the FULL dataset, the 95\% uncertainty-band widths at $x_{50}$ were 0.32 and 0.88 for the model-refitting and fixed-model analyses, respectively, under multiplicative noise, and 0.21 and 0.64 under additive noise. In the TRIM dataset, the corresponding widths were 0.69 and 1.00 under multiplicative noise and 0.98 and 0.97 under additive noise.

For both noise types, the uncertainty-band widths at $x_{50}$ were larger in the TRIM dataset than in the FULL dataset. The fixed-model band was wider than the model-refitting band in three of the four dataset--noise combinations. The exception was the TRIM dataset under additive noise, for which the two widths were similar (0.98 for model refitting and 0.97 for the fixed model). Within the fixed-model analysis, multiplicative noise produced wider bands than additive noise in both datasets.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\figWfull]{Figures/NOISEBAYES.pdf}
\caption{Measurement-noise propagation through the constrained Bayesian model for the FULL (A--B) and TRIM (C--D) datasets. Panels A and C show multiplicative Gaussian perturbations ($\sigma_{\mathrm{mult}}=0.129$), and Panels B and D show additive Gaussian perturbations ($\sigma_{\mathrm{add}}=0.144$). The black curve represents the model fitted to the unperturbed data. Dark and light shaded regions represent the 95\% uncertainty bands from the model-refitting and fixed-model analyses, respectively. The vertical dashed line indicates the clean-data midpoint $x_{50}$. Observed NC and AE data are shown at probabilities 0 and 1, respectively.}
\label{fig:noise_cb}
\end{figure}

For the constrained Bayesian model (Figure~\ref{fig:noise_cb}), the clean-data midpoint was $x_{50}=1.79$ for the FULL dataset and $x_{50}=1.75$ for the TRIM dataset. In the FULL dataset, the 95\% uncertainty-band widths at $x_{50}$ were 0.33 and 0.76 for the model-refitting and fixed-model analyses, respectively, under multiplicative noise, and 0.39 and 0.72 under additive noise. In the TRIM dataset, the corresponding widths were 0.69 and 0.99 under multiplicative noise and 1.00 and 0.94 under additive noise.

For both noise types and both analyses, the uncertainty-band widths at $x_{50}$ were larger in the TRIM dataset than in the FULL dataset. The fixed-model band was wider than the model-refitting band in three of the four dataset--noise combinations. The exception was the TRIM dataset under additive noise, where the model-refitting band was slightly wider (1.00 versus 0.94). Within the fixed-model analysis, multiplicative noise produced wider bands than additive noise in both datasets.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion}

The problem addressed in this study was how to estimate pneumonitis risk from a quantitative imaging biomarker when the available data are limited, and the biomarker itself is subject to measurement variability. This is particularly relevant in the present cohort, which included 58 patients but only five pneumonitis events. We therefore examined the biomarker--risk relationship using model formulations with different assumptions and evaluated uncertainty in model parameters, fitted risk functions, and derived risk-function characteristics. Additional analyses considered the influence of the identified normal-control observation, biomarker measurement noise, and parameter sensitivity, thereby linking uncertainty in model estimation and measurement to uncertainty in the resulting risk estimates.

The standard logistic models provided the simplest representation of the biomarker--risk relationship. Changing the predictor from \(x\) to \(\log x\) produced only small differences in \(x_{50}\) within each dataset, whereas the corresponding differences in \(s_{50}\) were larger. Thus, in the present data, the estimated location of the risk transition was less dependent on predictor formulation than its estimated sharpness. The monotone cubic logistic model was included as a more flexible extension of this relationship. Its fitted curves were broadly consistent with the increasing-risk pattern observed in standard logistic regression, while allowing additional variation in the shape of the transition. The likelihood profiles showed that the additional cubic parameters were not equally well determined by the available observations. Given the small number of pneumonitis events, this provides a useful illustration of the additional parameter estimation required when greater functional flexibility is introduced \citep{Riley2020}.

The constrained Bayesian model provided an additional formulation in which the NC and AE biomarker distributions were modeled separately, and conditional risk was derived using Bayes' theorem. In addition to the increased number of model parameters, monotonicity and tail conditions were imposed on the resulting risk function. Monotonicity —the assumption that increasing biomarker values produce decreasing estimated pneumonitis risk—is an assumption about the direction of the biomarker--risk relationship \citep{Dumbgen2024}. The tail conditions controlled the fitted risk function at very low and very high biomarker values. These conditions are model assumptions rather than properties established by the present data, particularly outside the observed biomarker range.

The model-evaluation results illustrate the balance between fit and model complexity among the models for which AIC and BIC were reported. The constrained Bayesian model achieved the highest Bernoulli log-likelihood in both datasets, whereas the log-transformed standard logistic model had the lowest AIC and BIC. Thus, the improvement in likelihood obtained with the constrained Bayesian formulation was accompanied by the penalty associated with its larger number of estimated parameters. These results should not be interpreted as identifying a single preferred model, because the formulations also differ in their assumptions and intended roles within the analysis. Rather, the comparison provides the basis for examining whether differences in model formulation and parameterization are subsequently reflected in uncertainty in the fitted risk functions and their derived characteristics.
Measurement noise provided an additional component of uncertainty beyond that associated with model estimation. In the present analysis, additive and multiplicative perturbations, with magnitudes determined from the preceding characterization of noise in the \(^{18}\)F-FDG PET/CT data, were propagated through the fitted models. For both the logistic regression and constrained Bayesian models, the uncertainty-band widths evaluated at the clean-data \(x_{50}\) were larger in the TRIM dataset than in the FULL dataset. In several TRIM analyses, the resulting bands covered a large proportion of the probability range. The fixed-model bands were wider than the corresponding model-refitting bands in three of the four datasets--noise combinations for each model. Within the fixed-model analysis, multiplicative perturbations also produced wider bands than additive perturbations in both datasets. These comparisons describe variability in predicted probability at the fixed \(x_{50}\) obtained from the unperturbed data and should therefore be distinguished from uncertainty in the estimated value of \(x_{50}\) itself. The effect of biomarker perturbation on predicted probability varied across the fitted risk function. Larger vertical changes in predicted probability were observed in regions where the fitted risk function changed more rapidly with the biomarker. This pattern was observed under both the logistic regression and constrained Bayesian formulations. The measurement-noise analysis, therefore, characterized a component of uncertainty distinct from that associated with model estimation. Together with the FULL--TRIM comparison, these results show that model-estimation uncertainty and uncertainty arising from perturbation of the measured biomarker did not necessarily exhibit the same pattern in the present analysis.

Uncertainty quantification provided a more informative distinction between uncertainty in model parameters and uncertainty in the fitted risk relationship. For standard logistic regression in the FULL dataset, the nonparametric and parametric bootstrap procedures produced extremely wide intervals for the regression coefficients, whereas the corresponding intervals for \(x_{50}\) remained on a substantially narrower biomarker range. The nonparametric bootstrap distribution was also multimodal, whereas the local Gaussian approximation captured the dominant mode near the maximum-likelihood estimate. Thus, in this small-event setting, resampling exposed parameter configurations not represented by the local approximation, but this large parameter variability did not propagate proportionally to the uncertainty in the estimated transition location. Parameter confidence intervals alone would therefore provide an incomplete characterization of epistemic uncertainty in the resulting risk assessment.

The constrained Bayesian model reinforced this distinction while also showing that uncertainty in a risk-function characteristic can depend strongly on the uncertainty-quantification procedure. Several model parameters had wide intervals under the Wald, NPBS, and PBS approaches. Although the point estimate of \(x_{50}\) was unchanged across these procedures within each dataset, its interval width was method-dependent, most notably for the PBS analysis of the TRIM dataset. At the level of the fitted risk function, the confidence bands were widest where the predicted probability changed most rapidly and were generally narrower for the TRIM than the FULL dataset. These results show that epistemic uncertainty is not represented by a single model-level quantity: parameter intervals characterize uncertainty in model estimation, confidence bands characterize its expression across the biomarker--risk function, and intervals for \(x_{50}\) characterize uncertainty specifically in the estimated transition location.

The aleatoric uncertainty considered in this study arose from measurement variability of the imaging biomarker and exhibited different FULL--TRIM behavior. For both the standard logistic and constrained Bayesian models, noise-induced uncertainty-band widths evaluated at the clean-data \(x_{50}\) were larger in the TRIM dataset than in the FULL dataset for both additive and multiplicative perturbations. This contrasts with the narrower bootstrap coefficient intervals observed after trimming for standard logistic regression and with the generally narrower model-based confidence bands for the constrained Bayesian model. Removal of the identified observation could therefore reduce some manifestations of epistemic uncertainty while increasing the propagation of biomarker measurement variability into predicted risk at the midpoint. The effect of an influential observation on model-estimation uncertainty should consequently not be assumed to predict its effect on aleatoric uncertainty in the resulting risk estimates.

The noise-propagation results also show that the magnitude of aleatoric uncertainty depended on both the perturbation model and whether the fitted model was held fixed or re-estimated. Fixed-model bands were wider than model-refitting bands in three of the four datasets--noise combinations for each model, and multiplicative perturbations produced wider bands than additive perturbations in both datasets within the fixed-model analysis. Several TRIM bands approached the full probability range at \(x_{50}\), demonstrating that the measured biomarker can introduce substantial variability in predicted risk even when uncertainty in the estimated location of \(x_{50}\) is considered separately. These bands quantify vertical variability in predicted probability at the fixed clean-data midpoint and should therefore not be interpreted as confidence intervals for \(x_{50}\).

Sensitivity analysis complements uncertainty analysis by identifying how strongly the characteristics of the derived risk function respond locally to changes in individual model parameters. This distinction is important because a parameter may be estimated with substantial uncertainty without having an equally large influence on \(x_{50}\) or \(s_{50}\), and conversely, a risk-function characteristic may be highly sensitive to a parameter even when parameter uncertainty is considered separately. Sensitivity and uncertainty, therefore, describe different properties of the fitted risk model and should not be interpreted interchangeably. For standard logistic regression, the parameter dependence of \(x_{50}\) was relatively balanced: within each formulation, the intercept and slope had equal absolute elasticities. In contrast, the dependence of \(s_{50}\) was strongly affected by predictor formulation. Under the RAW formulation, \(s_{50}\) was locally insensitive to the intercept and depended on the slope, whereas under the LOG formulation, both parameters contributed. The larger elasticities observed for \(s_{50}\) are consistent with the greater dependence of transition sharpness on model formulation observed in the fitted models. Thus, even when the RAW and LOG formulations locate the risk transition at similar biomarker values, the estimated sharpness of that transition is more dependent on the fitted parameterization. The constrained Bayesian model showed a more differentiated dependence on parameters. The prevalence parameter \(\omega\) had the smallest local influence on both \(x_{50}\) and \(s_{50}\), whereas the Gamma-distribution parameters \(k\) and \(\vartheta\) produced the largest absolute elasticities, with \(a\), \(b\), and \(s\) showing intermediate values. The absolute elasticities were also larger in the FULL than in the TRIM dataset for the parameters examined, with the difference particularly pronounced for \(s_{50}\). These results identify \(k\) and \(\vartheta\) as the parameters to which the derived risk-function characteristics were most locally responsive in the fitted Bayesian model. They also show that removal of the identified observation altered not only model estimates and uncertainty, but also the local dependence of the risk-function characteristics on the model parameters.
The sensitivity and uncertainty analyses provide complementary information about parameter behavior. Parameter intervals describe uncertainty in parameter estimation, whereas elasticity describes the local relative response of a risk-function characteristic to a proportional change in a parameter. The present results, therefore, distinguish parameters that are uncertain from those to which \(x_{50}\) and \(s_{50}\) are locally sensitive. More broadly, the combined analyses show that uncertainty in the fitted risk relationship cannot be inferred solely from model fit or parameter uncertainty: parameter uncertainty, risk-function uncertainty, noise-induced uncertainty, and local parameter sensitivity characterize distinct aspects of the resulting risk assessment.

Several limitations constrain the interpretation of these findings. First, the cohort was small, including only five patients with pneumonitis, which limited the precision of model estimates and the conclusions that could be drawn from comparisons among model formulations. The FULL--TRIM analysis further demonstrated the influence that a single observation can have in this setting. Second, the analysis was restricted to a single \(^{18}\)F-FDG PET/CT-derived biomarker, lung \(\SUV_{95\%}\), and therefore does not establish whether the findings extend to other SUV-based summaries, spatially resolved measures, or multi-organ biomarkers. Third, inference for the constrained Bayesian model was centered on the maximum a posteriori (MAP) estimate and used a local Laplace approximation rather than sampling from the full posterior distribution. The corresponding Wald intervals are therefore local approximations and should not be interpreted as a full posterior characterization; the NPBS and PBS provided complementary resampling-based assessments of uncertainty. Fourth, the aleatoric uncertainty considered here was limited to additive and multiplicative Gaussian perturbations of the derived subject-level biomarker. Although the preceding characterization informed the perturbation magnitudes of \(^{18}\)F-FDG PET/CT noise, other sources of imaging variability, including scanner-specific effects, reconstruction, segmentation, and patient motion, were not modeled separately.

Future work should evaluate the framework in larger, independent cohorts to assess predictive performance, calibration, and generalizability, and to determine whether the uncertainty behavior observed here persists with more outcome events. The framework could also be extended to alternative PET-derived biomarkers, longitudinal measurements, additional clinical variables, and multi-organ imaging features. Imaging-related uncertainty could be characterized further by evaluating segmentation and region-of-interest definitions and, where appropriate, spatially resolved lung measurements. Full posterior inference would also extend the MAP-based constrained Bayesian analysis by directly characterizing posterior uncertainty in the model parameters and derived risk-function characteristics.

The present study should therefore be interpreted as a methodological assessment of uncertainty-aware risk modeling rather than validation of a predictive biomarker for clinical use. Independent and prospective evaluation, followed by assessment of clinical utility, would be required before the resulting risk models could be considered for clinical decision-making.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}

As quantitative imaging biomarkers move from population-level associations to patient-specific risk estimation, uncertainty in the resulting estimates becomes an important component of model assessment, particularly when clinical data are limited and the biomarker is subject to measurement variability. In this study, we developed an uncertainty-aware quantitative risk assessment framework for modeling immune checkpoint inhibitor-related pneumonitis risk using \(^{18}\)F-FDG PET/CT-derived lung biomarkers.

The analysis showed that uncertainty in the fitted risk relationship cannot be inferred solely from model fit or parameter uncertainty. Large uncertainty in individual parameters did not propagate proportionally to uncertainty in the estimated transition location, and removal of the identified normal-control observation affected model-estimation uncertainty and noise-induced uncertainty differently. Measurement-noise propagation further showed that variability in the biomarker can produce substantial variability in predicted risk. At the same time, sensitivity analysis identified which model parameters had the greatest local influence on the location and sharpness of the risk transition. These findings demonstrate the value of assessing uncertainty at the levels of the model parameters, the fitted risk function, and derived characteristics of the risk function, while separately considering uncertainty introduced by biomarker measurement.

The present study is methodological and does not establish a clinically validated prediction model. The small number of pneumonitis events limits the conclusions that can be drawn from the fitted models. Larger independent cohorts are required to evaluate predictive performance, calibration, and generalizability and to characterize imaging-related uncertainty more comprehensively before the resulting risk models can be considered for clinical application.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ack 
The authors acknowledge the financial support from the Slovenian Research and Innovation Agency (research core funding No. P1-0389).

During the preparation of this manuscript, the authors used ChatGPT (OpenAI) to assist with language, clarity, and organization of the text. All scientific content, analyses, interpretations, and conclusions were reviewed and approved by the authors, who take full responsibility for the content of the manuscript.

\appendix
\renewcommand{\thesection}{Appendix \Alph{section}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Logistic regression models: formulation, fitting, and uncertainty quantification}
\label{sec:app_logistic}

This section summarizes the likelihood formulation, constrained extensions, and asymptotic approximations for the logistic regression models used in the main text.

% -------------------------------------------------------------------------
\subsection{Standard logistic regression}

\paragraph{Model.}
We consider observations of a positive predictor \(X\in\mathbb{R}_+\) and a binary response \(Y\in\{0,1\}\):
%

 \[
    \mathcal{D}
    =
    \{(x_i,y_i)\in\mathbb{R}_+\times\{0,1\}\}_{i=1}^n.
\]
%
To allow alternative parameterizations of the predictor, we define the transformation \(t_i=g(x_i)\), with \(g(x)=x\) (RAW) or \(g(x)=\log x\) (LOG), yielding.
%
\[
    \tilde{\mathcal{D}}
    =
    \{(t_i,y_i)\}_{i=1}^n.
\]
%
The conditional probability is modeled as
%
\[
  p(t;\bTheta)
  =
  \Prob(Y=1\mid X=x)
  =
  \expit\left(\eta_\text{LR}(t;\bTheta)\right),
\]
%
with a linear decision function
%
\[
  \eta_\text{LR}(t;\bTheta)
  =
  \bTheta^{\mathsf T}\boldsymbol{\phi}(t),
  \qquad
  \boldsymbol{\phi}(t)
  =
  (1,t)^{\mathsf T}.
\]

\paragraph{Likelihood.}
The Bernoulli log-likelihood is
%
\begin{equation}
    \ell_\text{LR}(\bTheta\mid\tilde{\mathcal{D}}) =
    \sum_{i=1}^n
    \Big[
      y_i\log p(t_i;\bTheta)
      +
      (1-y_i)\log\!\big(1-p(t_i;\bTheta)\big)
    \Big],
    \label{eq:bernoulli_llf}
\end{equation}
%
%with score
%
%\[
%  \nabla_{\bTheta}\ell_\text{LR}
%  =
%  \sum_{i=1}^n
%  (y_i-p(t_i;\bTheta))\boldsymbol{\phi}(t_i),
%\]
%
and its Hessian of the negative log-likelihood function is
%
\begin{equation}
  \mathbf{H}_\text{LR}(\bTheta) 
  = -\nabla_{\bTheta} \otimes \nabla_{\bTheta} \ell_\text{LR}(\bTheta)
  = \sum_{i=1}^n
  p(t_i;\bTheta)\big(1-p(t_i;\bTheta)\big)
  \boldsymbol{\phi}(t_i)\boldsymbol{\phi}(t_i)^{\mathsf T}.
  \label{eq:logit_hess}
\end{equation}

\paragraph{Asymptotic approximation.} The maximum likelihood estimator (MLE) is
%
\[
  \hat{\bTheta} =
  \argmax_{\bTheta\in\R^2}
  \ell_\text{LR}(\bTheta\mid\tilde{\mathcal{D}}).
\]
%
Under standard regularity conditions \citep{CasellaBerger2002},
%
\begin{equation}
  \bphi
  \mathrel{\dot{\sim}}
  \mathcal{N}
  \left(
    \hat{\bphi},
    \hat{\bSigma}_{\bphi}
  \right)
  \label{eq:bayes_phi_asymp}
\end{equation}
%
with a covariance matrix estimated as
%
\[
  \hat{\bSigma}_\text{LR} \approx \mathbf{H}_\text{LR}(\hat{\bTheta})^{-1}.
\]
%

Here, \(\mathbf{H}_\text{LR}(\hat{\bTheta})\) serves as an estimate of the Fisher information matrix, whose inverse provides the asymptotic covariance of the parameter estimates.

% -------------------------------------------------------------------------
\subsection{Monotone polynomial logistic models}

This model provides an intermediate level of flexibility between the standard linear logistic model and more structured constrained Bayesian models. To increase flexibility while preserving monotonicity, we replace the linear predictor by a polynomial in \(t\).

\paragraph{General formulation.}
We define
%
\begin{equation}
  p(t;\bTheta) =
  \expit\big(\eta_\text{poly LR}(t;\bTheta)\big),
  \label{eq:logit_p_gen}
\end{equation}
%
with
%
\begin{equation}
  \eta_\text{poly LR}(t;\bTheta) =
  \sum_{i=0}^p \theta_i t^i,
  \label{eq:logit_eta_gen}
\end{equation}
%
where \(\bTheta=(\theta_i)_{i=0}^p\). Since \(\expit\) is strictly increasing, monotonicity of \(p(t;\bTheta)\) is equivalent to requiring \(\eta_\text{poly LR}(t;\bTheta)\) to be non-decreasing. Therefore,
%
\[
  \eta_\text{poly LR}'(t;\bTheta) =
  \sum_{i=1}^p i\,\theta_i t^{i-1}
  \ge 0,
  \qquad
  \forall t\in\mathbb{R}.
\]
%
This defines the feasible region
%
\begin{equation}
  \Omega =
  \left\{
    \bTheta:
    \eta'(t;\bTheta)\ge 0
    \quad \forall t\in\mathbb{R}
  \right\}.
\end{equation}
%
The corresponding likelihood function is obtained by substituting Eq.~\eqref{eq:logit_p_gen} into Eq.~\eqref{eq:bernoulli_llf}, yielding \(\ell_\text{poly LR}(\bTheta\mid\tilde{\mathcal{D}})\). The estimator is then
%
\[
  \hat{\bTheta} =
  \arg\max_{\bTheta\in\Omega}
  \ell_\text{poly LR}(\bTheta\mid\tilde{\mathcal{D}}).
\]

\paragraph{Cubic logistic regression.}
For \(p=3\), we introduce the reparameterization
%
\begin{equation}
  \bTheta(\balpha)
  =
  \left(
    \alpha_0,\;
    \alpha_1^2+\alpha_3^2,\;
    \alpha_2\alpha_3,\;
    \frac{\alpha_2^2}{3}
  \right),
  \label{eq:cubic_transf}
\end{equation}
%
yielding
%
\begin{equation}
 \eta_\text{cubic LR}(t;\bTheta) =
  \alpha_0 + (\alpha_1^2+\alpha_3^2)t + (\alpha_2\alpha_3)t^2 +
  \frac{\alpha_2^2}{3}t^3,
  \label{eq:logit_cubic_odds}
\end{equation}
%
with derivative
%
\[
  \eta_\text{cubic LR}'(t;\bTheta) =
  (\alpha_3+\alpha_2 t)^2+\alpha_1^2
  \ge 0,
\]
%
for all \(\balpha\in\R^4\), ensuring global monotonicity by construction.

\paragraph{Estimation and asymptotics.}
The likelihood can be written as
%
\[
  \ell_{\text{cubic LR},\alpha}(\balpha) =
  \ell_\text{poly LR}\big(\bTheta(\balpha)\mid\tilde{\mathcal{D}}\big),
\]
%
with estimator
%
\[
  \hat{\balpha} =
  \arg\max_{\balpha\in\R^4}
  \ell_{\text{cubic LR},\alpha}(\balpha).
\]
%
This induces optimal parameters in original parametrization \( \hat{\bTheta} = \bTheta(\hat{\balpha})\). Under standard regularity conditions \citep{CasellaBerger2002},
%
\[
  \balpha
  \mathrel{\dot{\sim}}
  \mathcal{N}\!\left(
    \hat{\balpha},
    \bSigma_{\text{cubic LR},\balpha}
  \right),
\]
%
where the covariance is obtained from the inverse Hessian of the log-likelihood. The covariance of \(\bTheta\) is then obtained via the Wald approximation. Let
%
\[
  \mathbf{J}_\Theta(\balpha)
  =
  \frac{\partial \bTheta(\balpha)}
  {\partial\balpha^\top},
\]
%
denote the Jacobian of Eq.~\eqref{eq:cubic_transf}. Then
%
\[
  \bSigma_{\text{cubic LR},\bTheta}
  =
  \mathbf{J}_\Theta(\hat{\balpha})
  \bSigma_{\text{cubic LR},\balpha}
  \mathbf{J}_\Theta(\hat{\balpha})^\top.
\]

This construction guarantees global monotonicity of the fitted probability function while providing greater flexibility than the standard linear logistic model.

% -------------------------------------------------------------------------
\subsection{Confidence intervals and bands estimation for logistic regression}
\label{sec:app_ci_logistic}

This subsection provides the mathematical details for the uncertainty intervals of the logistic regression models introduced in the main text. We consider four complementary approaches: wald approximation, Monte-Carlo simulation propagation, nonparametric bootstrap, and parametric bootstrap.

\paragraph{Wald approximation.}

For the Wald approximation, uncertainty was propagated using the asymptotic normal approximation of the maximum-likelihood estimator (MLE) \citep{CasellaBerger2002}. Under standard regularity conditions,
%
\[
  \bTheta
  \mathrel{\dot{\sim}}
  \mathcal{N}
  \left(
    \hat{\bTheta},
    \hat{\bSigma}_{\mathrm{LR}}
  \right),
\]
%
where \(\hat{\bSigma}_{\mathrm{LR}}\) denotes the covariance matrix estimated from the inverse observed Fisher information.

For the fitted log-odds
%
\[
  \eta(t;\bTheta)=\bTheta^\top\boldsymbol{\phi}(t),
\]
%
where \(t=g(x)\) denotes the transformed predictor, the propagated variance at predictor value \(t\) is
%
\begin{equation}
  \mathrm{Var}\!\big[\eta(t;\hat{\bTheta})\big]
  =
  \boldsymbol{\phi}(t)^\top
  \hat{\bSigma}_{\mathrm{LR}}
  \boldsymbol{\phi}(t).
  \label{eq:eta_var_logistic}
\end{equation}
%
This yields two-sided confidence intervals with confidence level \(1-\alpha\) (significance level \(\alpha\)) for the log-odds evaluated at predictor value \(t\):
%
\begin{equation} 
  \eta(t;\hat{\bTheta})
  \pm
  z_{1-\alpha/2}
  \sqrt{
  \boldsymbol{\phi}(t)^\top
  \hat{\bSigma}_{\mathrm{LR}}
  \boldsymbol{\phi}(t)
  }.
\label{eq:eta_ci_logistic}
\end{equation}
%
where $z_p$ is the $p$th quantile of the normal distribution. The resulting interval bounds are transformed to the probability scale through the logistic link:
%
\[
  p(t;\bTheta)=\expit(\eta(t;\bTheta)).
\]

For the midpoint \(x_{50}\), uncertainty propagation is obtained via the Wald approximation:
%
\begin{equation}
  \mathrm{Var}(x_{50})
  \approx
  \left(
  \frac{\partial x_{50}(\hat{\bTheta})}{\partial \bTheta}
  \right)^\top
  \hat{\bSigma}_{\mathrm{LR}}
  \left(
  \frac{\partial x_{50}(\hat{\bTheta})}{\partial \bTheta}
  \right).
  \label{eq:x50_delta_var}
\end{equation}
%
Here,
%
\(
  \partial x_{50}(\hat{\bTheta})/\partial \bTheta
\)
%
denotes the gradient vector of the scalar function \(x_{50}(\bTheta)\) with respect to the parameter vector \(\bTheta\). Gradients are evaluated numerically using finite differences.

\paragraph{Monte-Carlo simulation propagation.}

To reduce reliance on the first-order local linearization inherent in the Wald approximation, parameter uncertainty was also propagated by direct sampling from the asymptotic distribution:
%
\begin{equation}
    \bTheta^{(j)} \sim
    \mathcal{N}(\hat{\bTheta},\hat{\bSigma}_{\mathrm{LR}}),
    \qquad
    j=1,\ldots,N.
\label{eq:normal_draws_logistic}
\end{equation}
%
Here, \(j\) denotes the index of the sampled parameter vector. For each draw, the corresponding fitted curve
%
\[
  p^{(j)}(t)=p(t;\bTheta^{(j)})
\]
%
and midpoint
%
\[
  x_{50}^{(j)}
\]
%
were recomputed. Confidence intervals were then obtained from empirical quantiles of the simulated distributions. Unlike the Wald approximation, this approach propagates uncertainty through repeated sampling rather than local linearization.

\paragraph{Nonparametric bootstrap.}

Bootstrap dataset 

%
\[
  \mathcal D^{(j)}, \qquad j=1,\ldots,N,
\]
%
were generated by stratified resampling with replacement within AE and NC groups, preserving class counts, where \(N\) denotes the total number of bootstrap replicates and $j$ is its index. Each bootstrap dataset was refitted to obtain parameter and midpoint estimates:
%
\[
  \hat{\bTheta}^{(j)}, \qquad x_{50}^{(j)}.
\]
%
Confidence intervals were constructed from empirical quantiles of the resulting bootstrap distributions. This procedure captures sampling variability directly from the observed cohort without assuming a parametric data-generating model.

\paragraph{Parametric bootstrap.} For the parametric bootstrap, predictor values \(t_i = g(x_i)\) were kept fixed, while outcomes were simulated from the fitted Bernoulli model:
%
\begin{equation}
  Y_i^{(j)} \sim \mathrm{Bernoulli}\!\left(p(t_i;\hat{\bTheta})\right),
  \qquad
  i=1,\ldots,n,\quad j=1,\ldots,N,
  \label{eq:param_boot_logistic}
\end{equation}
%
where \(n\) denotes the number of patients and \(N\) denotes the total number of bootstrap replicates. Each simulated dataset was refitted to obtain
%
\[
  \hat{\bTheta}^{(j)}, \qquad x_{50}^{(j)}.
\]
%
Confidence intervals were again computed from empirical quantiles across bootstrap replicates. This procedure quantifies uncertainty under the fitted logistic model structure.

\paragraph{Primary inferential target.} Although confidence bands for the fitted risk function were computed, the principal inferential quantity was the midpoint \(x_{50}\). Accordingly, confidence intervals for \(x_{50}\) were obtained directly from its propagated or resampled distribution, rather than inferred indirectly from confidence bands on the model value.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Constrained Bayesian model: formulation, fitting, and uncertainty quantification}
\label{app:bayesian_development}

This appendix presents the mathematical development of the constrained Bayesian model introduced in Section~\ref{sec:bayesian_model}. We first define the Bayesian risk model, then discuss the selection of group-conditional distribution families, and finally present the structural constraints and optimization procedure used to determine the optimal model parameters.

% -------------------------------------------------------------------------
\subsection{Model definition}
\label{app:conditional-bayes-model}

We consider a random pair \((X,Y)\), where \(X \ge 0\) denotes the biomarker value and \(Y \in \{\NCgrp,\AEgrp\}\) indicates the patient group: normal control \((\NCgrp)\) or adverse event \((\AEgrp)\). In the calculations, we use the numerical encoding \(y_i=0\) for \(\NCgrp\) and \(y_i=1\) for \(\AEgrp\). A sample of random pairs,
%
\[
\mathcal D
=
\{(x_i,y_i)\}_{i=1}^{n},
\]
%
constitutes our dataset.

Following the discussion in Section~\ref{sec:bayesian_model}, the conditional probability of an adverse event,
%
\[
  P(Y=\AEgrp\mid X=x),
\]
%
is referred to as the risk. By Bayes' rule, it is given by
%
\begin{equation}    
  p(x;\bTheta)=
  \frac{\omega\,p_{X\mid Y}(x\mid\AEgrp)}
  {\omega\,p_{X\mid Y}(x\mid\AEgrp)+(1-\omega)\,p_{X\mid Y}(x\mid\NCgrp)},
  \label{eq:risk_fun}
\end{equation}
%
where \(\omega=P(Y=\AEgrp)\) denotes the marginal prevalence of the \(\AEgrp\) group. This expression can be rewritten into logistic form:
%
\[
  p(x;\bTheta)=\frac{1}{1+\psi(x;\bTheta)},
\]
%
with
%
\[
  \psi(x;\bTheta)=
  \frac{1-\omega}{\omega}
  \frac{p_{X\mid Y}(x\mid\NCgrp)}
       {p_{X\mid Y}(x\mid\AEgrp)}.
\]
%
Equivalently,
%
\[
  p(x;\bTheta)=\expit\!\big(\eta_{\mathrm{CB}}(x;\bTheta)\big),
\]
%
where
%
\[
  \eta_{\mathrm{CB}}(x;\bTheta)
  =
  -\log\psi(x;\bTheta).
\]
%

Thus, \(\eta_{\mathrm{CB}}(x;\bTheta)\) plays the role of the decision function in the constrained Bayesian model.

The precise functional forms of these distributions are determined in the next subsection through distribution-family selection and structural considerations.

% -------------------------------------------------------------------------
\subsection{Selection of distribution families}
\label{sec:app_families}

The choice of class-conditional biomarker distributions plays a central role in shaping the behavior of the Bayesian risk function. Model behavior depends not only on the adequacy of marginal fits within each group, but also on the structure of the density ratio appearing in Bayes' rule. Given that the biomarker is strictly positive and right-skewed, we examined several positive-support candidate distributions, including Gamma, LogNormal, and Beta-prime. Distribution selection was guided by three main criteria:
%
\begin{itemize}
\item[i)] within-group fit,
\item[ii)] stability and monotonicity of the induced risk function,
\item[iii)] consistency of asymptotic and tail behavior.
\end{itemize}
%
To obtain a conditional probability that increases from 0 to 1 as the biomarker value increases, one may use more flexible distributions defined on \((0,\infty)\). Such distributions often involve additional parameters, which can be difficult to justify in small datasets because the available data may not reliably estimate all of them. We first illustrate these issues using Gaussian-type models and then consider positive-support alternatives.

% -------------------------------------------------------------------------
\subsubsection{Gaussian-type distributions}
\label{sec:app_families_gauss}

We first examined Gaussian-type class-conditional distributions as a natural baseline choice for modeling biomarker distributions in the Bayesian framework. This model is referred to as the Gaussian Bayesian (GB) model.

Let the group-conditional biomarker distributions be Gaussian:
%
\begin{equation}
  p_{X\mid Y}(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\{\AEgrp,\NCgrp\},
  \label{eq:gauss_grp_cond}
\end{equation}
%
where \(\mu_y\) and \(\sigma_y\) denote the group-specific mean and standard deviation. Using Bayes' theorem, the conditional probability can be written in logistic form:
%
\[
  \Prob(Y=\AEgrp\mid X=x)
  =
  \expit(\eta_\text{GB}(x;\bTheta)),
\]
%
with a decision function
%
\begin{equation}
  \eta_\text{GB}(x;\bTheta)
  =
  \ln\frac{\omega}{1-\omega}
  +
  \ln\frac{\sigma_{\NCgrp}}{\sigma_{\AEgrp}}
  -
  \frac{(x-\mu_{\AEgrp})^2}{2\sigma_{\AEgrp}^2}
  +
  \frac{(x-\mu_{\NCgrp})^2}{2\sigma_{\NCgrp}^2},
  \label{eq:gauss_eta}
\end{equation}
%
where
%
\[
 \bTheta=
  (\omega,\mu_{\NCgrp},\sigma_{\NCgrp},\mu_{\AEgrp},\sigma_{\AEgrp}).
\]
%
Importantly, the decision function is quadratic in \(x\):
%
\[
 \eta_\text{GB}(x;\bTheta) = a_0+a_1x+a_2x^2.
\]
%
Thus, although the parameter vector contains five parameters, the resulting decision function depends effectively on only three coefficients, leading to a partially non-identifiable parameterization.

We distinguish two cases:

\begin{itemize}

\item[i)]
If \(\sigma_{\NCgrp}\neq\sigma_{\AEgrp}\), the quadratic term dominates the asymptotic behavior:
%
\begin{equation}
\lim_{x\to\pm\infty}\Prob(Y=\AEgrp\mid X=x) = \left\{
\begin{array}{lll}
1&:& \sigma_{\NCgrp}<\sigma_{\AEgrp},\\
0&:& \sigma_{\NCgrp}>\sigma_{\AEgrp}.
\end{array}\right.
\end{equation}
%
Thus, the Gaussian specification does not generally yield a globally stable monotone S-shaped risk function.

\item[ii)]
If \(\sigma_{\NCgrp}=\sigma_{\AEgrp}=\sigma\), the model reduces to a linear log-odds form:
%
\[
\Prob(Y=\AEgrp\mid X=x) =
\frac{1}{1+c\exp(-b(x-a))},
\]
%
with
%
\[
a=\frac{\mu_{\NCgrp}+\mu_{\AEgrp}}{2},
\qquad
b=\frac{\mu_{\AEgrp}-\mu_{\NCgrp}}{\sigma^2},
\qquad
c=\frac{1-\omega}{\omega}.
\]
%
Monotonicity requires \(b>0\).

\end{itemize}

Overall, Gaussian-type distributions may lead to unstable or non-monotone behavior, particularly in the tails. This motivates the use of alternative positive-support distributions.

% -------------------------------------------------------------------------
\subsubsection{Positive-support candidate families}
\label{sec:app_families_candidates}

Based on the criteria in the first paragraph of \ref{sec:app_families} above, we examined several positive-support combinations for the class-conditional densities, including LogNormal--LogNormal, Gamma--Gamma, LogNormal--Gamma, and BetaPrime--Gamma specifications.

The log-normal models may fit each group reasonably well, but the resulting log-density ratio is quadratic in \(\log x\) and follows a Gaussian distribution. This can lead to unstable extrapolation or non-monotonic risk estimates outside the observed biomarker range. More generally, using the same distributional family for both groups may give acceptable goodness-of-fit results, but still produce undesirable behavior in the tails of the induced risk function.

Some distributions on the positive real line allow the density ratio in Bayes' rule to remain monotone for \(x>0\), while also giving better control over tail
behavior and asymptotic limits. To illustrate how such distributions can produce
a monotonically increasing sigmoid-like risk function \(p(x;\bTheta)\) in
Eq.~\eqref{eq:risk_fun}, we briefly consider two representative examples:
%
\begin{enumerate}

\item Beta-prime distributions \citep{Johnson1995}: Let the group-conditional biomarker distributions follow Beta-prime (BP) distributions
%
\[
p_{X\mid Y}(x\mid y) =
\frac{
\left(\frac{x}{s_y}\right)^{\alpha_y-1}
\left(1+\frac{x}{s_y}\right)^{-(\alpha_y+\beta_y)}
}
{s_y\,\mathrm{B}(\alpha_y,\beta_y)},
\qquad
y\in\{\AEgrp,\NCgrp\},
\quad x>0.
\]

%
where \(\alpha_y>0\) and \(\beta_y>0\) are shape parameters, \(s_y>0\) is a scale parameter, and \(\mathrm{B}(\alpha_y,\beta_y)\) denotes the beta function. The group-conditional densities given by
%
\[
  X \mid (Y=y)
  \sim
  \mathrm{BetaPrime}(\alpha_y,\beta_y,s_y),
  \qquad
  y\in\{\AEgrp,\NCgrp\}.
\]
%
Under suitable parameter conditions, the induced risk function may exhibit monotone-increasing sigmoid-like behavior. For example, this occurs when%
\[
\alpha_{\AEgrp}\ge \alpha_{\NCgrp},
\qquad
\beta_{\NCgrp}\ge \beta_{\AEgrp},  
\qquad
(\alpha_{\AEgrp}+\beta_{\NCgrp})s_{\AEgrp}
\ge
(\alpha_{\NCgrp}+\beta_{\AEgrp})s_{\NCgrp}.
\]

\item Gamma distributions \citep{Johnson1994}.

Assume that the group-conditional biomarker distributions follow Gamma distributions:
%
\[
p_{X\mid Y}(x\mid y) =
\frac{
x^{k_y-1}e^{-x/\vartheta_y}
}{
\Gamma(k_y)\vartheta_y^{k_y}
},
\qquad
y\in\{\AEgrp,\NCgrp\},
\quad x>0.
\]
%
where \(k_y>0\) and \(\vartheta_y>0\) denote the shape and scale parameters, respectively, and \(\Gamma(k_y)\) denotes the gamma function.

Under suitable parameter relations, the induced risk function
%
\[
\Prob(Y=\AEgrp\mid X=x)
\]
%
may exhibit monotone increasing sigmoid-like behavior. For example, monotonicity may arise when
%
\[
k_{\NCgrp}<k_{\AEgrp},
\qquad
\vartheta_{\NCgrp}<\vartheta_{\AEgrp}.
\]

\end{enumerate}

These examples illustrate that several positive-support distribution families can generate increasing sigmoid-like risk profiles under appropriate parameter constraints. However, assigning independent distributional parameters to both groups rapidly increases model complexity. Even relatively simple specifications require multiple distributional parameters beyond the prevalence parameter. In small cohorts, such parameterizations may become weakly identifiable and statistically unstable.

Although several positive-support specifications were theoretically feasible, the BetaPrime--Gamma combination provided the most stable and interpretable behavior in our dataset. Empirically, the AE group exhibited a broader distribution with heavier upper tails, whereas the NC group was concentrated at lower biomarker values with a lighter right tail. The Beta-prime distribution naturally accommodates heavier right-tail behavior, while the Gamma distribution provides a parsimonious and stable representation of the lower-valued NC distribution. This combination also yielded a smooth monotone risk function with appropriate asymptotic behavior over the observed biomarker range.

Based on these considerations, the final model specification uses a Beta-prime distribution for the AE group and a Gamma distribution for the NC group:
%
\[
\begin{aligned}
X\mid(Y=\AEgrp)&\sim\operatorname{BetaPrime}(a,b,s),\\
X\mid(Y=\NCgrp)&\sim\operatorname{Gamma}(k,\vartheta).
\end{aligned}
\]
%
The corresponding group-conditional densities are
%
\begin{equation}
p_{X\mid Y}(x\mid\AEgrp)
=
\frac{
\left(\frac{x}{s}\right)^{a-1}
\left(1+\frac{x}{s}\right)^{-(a+b)}
}
{s\,\mathrm{B}(a,b)},
\qquad x>0,
\label{eq:bp_final_density}
\end{equation}
%
and
%
\begin{equation}
p_{X\mid Y}(x\mid\NCgrp)
=
\frac{
x^{k-1}e^{-x/\vartheta}
}{
\Gamma(k)\vartheta^k
},
\qquad x>0.
\label{eq:gamma_final_density}
\end{equation}

This choice accommodates the positive, right-skewed biomarker distribution, captures heavier upper tails in the AE group, and supports a monotone, asymptotically stable risk function.

Figure~\ref{fig:SuppDistributions} illustrates the empirical distributions together with the selected fitted parametric distributions.

\begin{figure}[!htbp]
\centering
\includegraphics[width=\linewidth]{Figures/Distributions.png}
\caption{Empirical biomarker distributions and fitted group-conditional densities. Histograms of lung \( \SUV_{95\%} \) for the NC and AE groups. The horizontal axis shows the biomarker value \(x=\SUV_{95\%}\), and the vertical axis shows the probability density function (PDF), representing normalized densities with unit area. Dashed curves indicate kernel density estimates, and solid curves show fitted parametric distributions. The final constrained Bayesian model uses a Gamma distribution for the NC group and a Beta-prime distribution for the AE group.}
\label{fig:SuppDistributions}
\end{figure}

% -------------------------------------------------------------------------
\subsubsection{Asymptotic behavior and monotonicity}
\label{sec:app_bayes_conditions}

Having selected Beta-prime and Gamma distributions for the \(\AEgrp\) and \(\NCgrp\) groups, respectively, we now derive the structural constraints required to ensure appropriate asymptotic behavior and monotonicity of the Bayesian risk function.

The model parameters are collected in the vector
%
\[
  \bTheta=(\omega,a,b,s,k,\vartheta),
\]
%
with natural parameter space \(\omega\in(0,1)\) and \(a,b,s,k,\vartheta>0\).

Since \(\expit\) and \(\logit\) are inverse bijections, it is sufficient to study the logit transformation of the risk function:
%
\[
  \logit\left(p(x;\bTheta)\right) =
  \logit(\omega)
  +
  \log p_{X\mid Y}(x\mid\AEgrp)
  -
  \log p_{X\mid Y}(x\mid\NCgrp).
\]
%
Substituting the class-conditional log-densities gives
%
\begin{align}
  \log p_{X\mid Y}(x\mid\AEgrp) &=
  (a-1)\log x
  -
  (a+b)\log\left(1+\frac{x}{s}\right)
  -
  \log\big(s\,\mathrm{B}(a,b)\big),
  \\
  \log p_{X\mid Y}(x\mid\NCgrp) &=
  (k-1)\log x
  -
  \frac{x}{\vartheta} 
  -
  \log\big(\Gamma(k)\vartheta^k\big),
\end{align}
%
which yields the discriminant term
%
\[
  \Delta(x;a,b,s,k,\vartheta)
  =
  (a-k)\log x
  -
  (a+b)\log\left(1+\frac{x}{s}\right)
  +
  \frac{x}{\vartheta}.
\]
%
Thus, the Bayesian decision function can be written as
%
\[
  \eta_{\mathrm{CB}}(x;\bTheta)
  =
  \logit(\omega)
  +
  \Delta(x;a,b,s,k,\vartheta),
\]
%
and the corresponding risk function is
%
\[
  p(x;\bTheta)
  =
  \expit\!\big(
  \eta_{\mathrm{CB}}(x;\bTheta)
  \big).
\]

We now derive conditions on the parameters that guarantee the required asymptotic behavior and monotonicity.

\begin{itemize}

\item[i)] To satisfy the asymptotic requirements
%
\[
  \lim_{x\to0^+}p(x;\bTheta)=0,
  \qquad
  \lim_{x\to\infty}p(x;\bTheta)=1,
\]
%
it is sufficient that
%
\[
  \lim_{x\to0^+}\Delta(x;\cdot)=-\infty,
  \qquad
  \lim_{x\to\infty}\Delta(x;\cdot)=\infty.
\]

These asymptotic conditions are satisfied when
%
\begin{equation}
  a>k,
  \qquad
  \vartheta>0.
  \label{eq:bayes_cond1}
\end{equation}
The second condition is already part of the positivity constraint for the Gamma scale parameter.

\item[ii)] Since the logistic function \(\expit\) is strictly increasing, the risk function \(p(x;\bTheta)\) is increasing whenever the log-density-ratio term is increasing. Therefore, we require
%
\[
  \Delta'(x;\cdot)>0,
  \qquad
  \text{for all }x>0.
\]

A sufficient global monotonicity condition is
%
\begin{equation}
  \vartheta\le\vartheta_{\max}(a,b,k,s),
  \qquad
  \vartheta_{\max}(a,b,k,s)
  =
  \frac{s}
  {\left(\sqrt{a+b}-\sqrt{a-k}\right)^2}.
  \label{eq:bayes_cond2}
\end{equation}

\end{itemize}

Under the constraints in Eqs.~\eqref{eq:bayes_cond1} and~\eqref{eq:bayes_cond2}, the Bayesian risk function \(p(x;\bTheta)\) is strictly increasing on \((0,\infty)\) and satisfies the desired endpoint limits:
%
\[
  \lim_{x\to0^+}p(x;\bTheta)=0,
  \qquad
  \lim_{x\to\infty}p(x;\bTheta)=1.
\]

These conditions define the admissible parameter region used in the constrained Bayesian model.
% -------------------------------------------------------------------------
\subsection{Objective function and parameter estimation} \label{sec:app_bayes_objfun}

We now describe how the optimal Bayesian parameters are determined within the admissible parameter region defined above.

For the constrained Bayesian (CB) risk model, the Bernoulli log-likelihood for the observed data is:
%
\[
  \ell_\text{CB}(\bTheta \mid \mathcal{D}) =
  \sum_{i=1}^n
  \left[
    y_i\log p(x_i;\bTheta)
    +
    (1-y_i)\log\{1-p(x_i;\bTheta)\}
  \right].
\]
%
where the sum runs over all observations in the dataset \(\mathcal{D}=\{(x_i,y_i)\}_{i=1}^n\). In a pure maximum likelihood approach, the parameter estimate would be obtained by maximizing \(\ell_\text{CB} (\bTheta \mid \mathcal{D})\). However, in practice, this optimization can be unstable because the likelihood may behave irregularly near the boundary values of the prevalence parameter \(\omega\). To stabilize estimation, we place a beta prior on \(\omega\):
%
\[
  \omega\sim\operatorname{Beta}(\alpha,\beta).
\]
%
The parameters \(\alpha\) and \(\beta\) are chosen so that the prior mean equals the empirical (sample) prevalence
%
\[
  \omega_{\mathrm{emp}} =\frac{1}{n}\sum_{i=1}^n y_i.
\]
%
Specifically, we set
%
\[
  \alpha = \tau\omega_{\mathrm{emp}},
  \qquad
  \beta = \tau(1-\omega_{\mathrm{emp}}),
\]
%
where \(\tau>0\) controls the concentration of the prior. With this choice,
%
\[
  \Expect[\omega]  = \omega_{\mathrm{emp}},
  \qquad
  \Var[\omega] = \frac{\omega_{\mathrm{emp}}(1-\omega_{\mathrm{emp}})}{\tau+1}.
\]
%
For \(\tau\), we do not need a precise estimate, but rather a rough estimate that allows for natural variation in \(\omega\) across samples. Based on the previous studies summarized in Table~\ref{tab:irae_clinical_trials}, we find that the reported values of \(\omega\) range from \(0\%\) to \(9\%\), with mean and sample standard deviation \(4.1\% \pm 3.1\%\). Our empirical estimate is \(\omega_{\mathrm{emp}} \doteq 9\%\). Based on these results, we assume that the standard deviation of \(\omega\) is approximately \(5\%\), which yields \(\tau \approx 30\). In the spirit of Bayesian estimation, we define the objective function as the sum of the Bernoulli log-likelihood and the log-density of the beta prior:
%
\begin{equation}
  L_\text{CB}(\bTheta\mid \mathcal{D}) = 
  \ell_\text{CB}(\bTheta\mid \mathcal{D}) + (\alpha-1)\log\omega + (\beta-1)\log(1-\omega).
  \label{eq:bayes_objfun}
\end{equation}
%
In Bayesian statistics, this objective function corresponds to the log-posterior density up to an additive constant \citep{CasellaBerger2002, Gelman2013}. The parameter estimate is therefore obtained by solving
%
\[
  \hat{\bTheta}
  =
  \argmax_{\bTheta\in\Omega}
  L_\text{CB}(\bTheta \mid \mathcal{D}),
\]
%
where \(\Omega\) denotes the admissible parameter region defined by Eqs.~\eqref{eq:bayes_cond1} and~\eqref{eq:bayes_cond2}. Because of this connection with Bayesian estimation, we refer to \(\hat{\bTheta}\) as the maximum a posteriori (MAP) estimate.


\begin{table}[htbp]
\centering
\caption{Percents of immune-related adverse events reported in selected melanoma immunotherapy clinical trials from references: [1] \citep{Hribernik2022}, [2]\citep{Robert2021}, [3] \citep{Robert2019}, [4] \citep{Larkin2018}, [5]\citep{Wolchok2022}, [6] \citep{Hodi2016}.}
\label{tab:irae_clinical_trials}

\begin{tabularx}{\textwidth}{@{}llXcccc@{}}
\toprule
Treatment & Colitis & Hypothyroidism & Hyperthyroidism & Pneumonitis & Reference\\
\midrule
H anti-PD-1 or anti-CTLA-4 & 10\% & 16\% & 16\% & 9\% & [1]\\
Pembrolizumab (anti-PD-1) & 2.3\% & 9\% & 3\% & 3\% &  [2] \\
Pembrolizumab (anti-PD-1) & 2.3\% & 9\% & 3\% & 3\% & [2] \\
Pembrolizumab (anti-PD-1) & 2.3\% & 9\% & 3\% & 3\% & [2] \\
Nivolumab (anti-PD-1) & 2\% & 7\% & 2\% & 3\% &  [3]\\
Nivolumab (anti-PD-1) & 3\% & 11\% & 4\% & 2\%  & [3] \\
Ipilimumab (anti-CTLA-4) & 11\% & 5\% & 1\% & 2\%  & [5]\\
Ipilimumab + Nivolumab & 13\% & 18\% & 12\% & 7\%  & [5]\\
Ipilimumab (anti-CTLA-4) & 6\% & 13\% & 0\% & 0\%  & [6]\\
Ipilimumab + Nivolumab & 18\% & 17\% & 2\% & 9\%  &  [6] \\
\bottomrule
\end{tabularx}
\end{table}

% -------------------------------------------------------------------------
\subsection{Unconstrained parameterization and transformations} 
\label{sec:app_bayes_transf}

Direct optimization under nonlinear parameter constraints is computationally challenging. To perform unconstrained optimization, we introduce an internal parameter vector
%
\[
  \bphi=
  (\phi_\omega,\phi_b,\phi_s,\phi_k,\phi_\delta,\phi_\vartheta)
  \in\R^6.
\]
%
The transformation
%
\[
  \mathbf{g}:\bphi\mapsto\bTheta
\]
%
maps the unconstrained optimization parameters to the constrained model parameter space and is defined by
%
\[
  \begin{split}
    \omega &= \expit(\phi_\omega), \qquad
    b = \softplus(\phi_b), \qquad
    s = \softplus(\phi_s), \\
    k &= \softplus(\phi_k), \qquad
    \delta = \softplus(\phi_\delta), \qquad
    a = k+\delta.
  \end{split}
\]
%
where \( \softplus(x)=\log(1+\exp(x))\). Since \(\delta>0\), the constraint
%
\[
a>k
\]
%
from Eq.~\eqref{eq:bayes_cond1} is automatically satisfied. We then set
%
\[
\vartheta =
  \vartheta_{\max}(a,b,k,s)\expit(\phi_\vartheta).
\]
%
This guarantees
%
\[
  0<\vartheta<\vartheta_{\max}(a,b,k,s),
\]
%
so the monotonicity constraint in Eq.~\eqref{eq:bayes_cond2} is satisfied automatically. The optimizer therefore operates over the unrestricted vector \(\bphi\in\R^6\), while the model is evaluated using the constrained parameter vector
%
\[
  \bTheta=\mathbf{g}(\bphi).
\]
%
The MAP estimate is obtained by optimizing the objective function \(L_\text{CB}\), Eq.~\eqref{eq:bayes_objfun}, in the internal parameter space:
%
\[
  \hat{\bphi} = \argmax_{\bphi\in\R^6} 
  L_\text{CB}\!\left(\mathbf{g}(\bphi)\mid\mathcal{D}\right),
\]
%
and mapping the resulting estimate back to the original parameter space: \( \hat{\bTheta} = \mathbf{g}(\hat{\bphi})\). This smooth reparameterization ensures that every parameter vector explored during optimization automatically satisfies the structural constraints required for valid probabilities, monotonicity, and stable tail behavior, while allowing efficient unrestricted optimization in Euclidean space.

%-------------------------------------------------------------------------------------------------------------------
\subsection{Confidence interval and band estimation for the constrained Bayesian model}
\label{sec:app_ci_bayes}

This subsection summarizes the construction of uncertainty intervals for the constrained Bayesian model introduced in the main text. In contrast to logistic regression, inference is based on the maximum a posteriori (MAP) estimate, and uncertainty is approximated locally using curvature of the log-posterior, together with resampling-based approaches.

\paragraph{Wald approximation.}

Let \(\bphi\in\mathbb{R}^p\) denote the unconstrained parameter vector used for optimization, and let
%
\[
  \bTheta=\mathbf g(\bphi)
\]
%
denote the corresponding constrained parameterization, where \(\mathbf g(\cdot)\) is the smooth parameter-space transformation introduced in the unrestricted reparameterization subsection. 

The MAP estimate of the optimization parameters is defined as
%
\[
  \hat{\bphi}
  =
  \arg\max_{\bphi}
  L_{\mathrm{CB}}(\mathbf g(\bphi)\mid \mathcal{D}),
  \qquad
  \hat{\bTheta}
  =
  \mathbf g(\hat{\bphi}),
\]
%
where \(\mathcal D\) denotes the observed dataset and \(L_{\mathrm{CB}}\) is the constrained Bayesian objective function defined in~\ref{sec:app_bayes_objfun}.

Under standard regularity conditions \citep{CasellaBerger2002}, the local quadratic approximation of the objective function  around \(\hat{\bphi}\) yields an approximate distribution of parameters
%
\[
  \bphi
  \mathrel{\dot{\sim}}
  \mathcal{N}
  \left(
    \hat{\bphi},
    \hat{\bSigma}_{\bphi}
  \right),
\]
%
with covariance matrix approximated by the inverse observed Hessian:
%
\[
  \hat{\bSigma}_{\bphi}
  \approx
  \mathbf{H}(\hat{\bphi})^{-1},
\]
%
where
%
\[
  \mathbf{H}(\hat{\bphi})
  =
  -\left.
  (\nabla_{\bphi}\otimes\nabla_{\bphi})
  L_{\mathrm{CB}}(\mathbf g(\bphi)\mid \mathcal{D})
  \right|_{\bphi=\hat{\bphi}}.
\]
%
denotes the Hessian matrix evaluated at \(\hat{\bphi}\). This local Gaussian approximation provides an analytical approximation of parameter uncertainty near the MAP estimate.

\paragraph{Transformation to model parameters.}

Uncertainty is propagated to the constrained parameters \(\bTheta\) through the transformation
%
\[
 \bTheta=\mathbf g(\bphi).
\]
%
Let
%
\[
  \mathbf{J}_\Theta(\bphi) =
  \frac{\partial \mathbf g(\bphi)}{\partial \bphi^\top}
\]
%
denote the corresponding Jacobian. The covariance matrix of \(\bTheta\) is then approximated as
%
\begin{equation}
  \hat{\bSigma}_{\bTheta} =
  \mathbf{J}_\Theta(\hat{\bphi})
  \hat{\bSigma}_{\bphi}
  \mathbf{J}_\Theta(\hat{\bphi})^\top.
  \label{eq:bayes_theta_cov}
\end{equation}

\paragraph{Propagation to model outputs.}

Let \(p(x;\bTheta)\) denote the fitted risk function. Its variance is approximated by
%
\begin{equation}
  \mathrm{Var}\!\big[p(x;\hat{\bTheta})\big]
  \approx
  \left(
    \frac{\partial p(x;\hat{\bTheta})}{\partial \bTheta}
  \right)^\top
  \hat{\bSigma}_{\bTheta}
  \left(
  \frac{\partial p(x;\hat{\bTheta})}{\partial \bTheta}
  \right).
  \label{eq:bayes_curve_var}
\end{equation}
%
with \(\partial p(x;\hat{\bTheta})/\partial \bTheta\) denoting the gradient vector of the scalar function \(p(x;\bTheta)\), evaluated at fixed \(x\), with respect to the parameter vector \(\bTheta\). Gradients are computed numerically using finite differences. Confidence bands with confidence level $1-\alpha$ are then obtained as
%
\[
  p(x;\hat{\bTheta}) \pm z_{1-\alpha/2} \sqrt{\mathrm{Var}[p(x;\hat{\bTheta})]},
\]
%
with truncation to the interval \([0,1]\), where $z_p$ is $p$th quantile of standard normal distribution.

The midpoint \(x_{50}\) is defined implicitly by \(p(x_{50};\bTheta)=0.5\). Its variance is approximated by
%
\begin{equation}
  \mathrm{Var}(x_{50})
  \approx
  \left(
    \frac{\partial x_{50}(\hat{\bTheta})}{\partial \bTheta} 
  \right)^\top
  \hat{\bSigma}_{\bTheta}
  \left(
    \frac{\partial x_{50}(\hat{\bTheta})}{\partial \bTheta}
  \right).
  \label{eq:bayes_x50_delta}
\end{equation}
%
where \(\partial x_{50}(\hat{\bTheta})/\partial \bTheta \) denotes the gradient vector of the scalar function \(x_{50}(\bTheta)\) with respect to the parameter vector \(\bTheta\). Gradients are computed numerically using finite differences.

\paragraph{Nonparametric bootstrap.}

Bootstrap datasets
%
\[
  \mathcal D^{(j)} =
  \left\{
    \big(x_i^{(j)},y_i^{(j)}\big)
  \right\}_{i=1}^{n},
  \qquad
  j=1,\ldots, N,
\]
%
where \(n\) denotes the number of patients and \(N\) denotes the total number of bootstrap replicates. The bootstrap datasets were generated by resampling the observed dataset \(\mathcal{D}\) with replacement. Each bootstrap dataset was refitted via MAP estimation to obtain
%
\[
  \hat{\bTheta}^{(j)}, \qquad x_{50}^{(j)}.
\]

Confidence intervals were constructed from empirical quantiles of the resulting bootstrap distributions. This procedure captures sampling variability directly from the observed cohort without imposing additional distributional assumptions.

\paragraph{Parametric bootstrap.}

For the parametric bootstrap, synthetic datasets were generated in the joint \((X,Y)\) space using the fitted generative model. Group labels were sampled as
%
\begin{equation}
  Y_i^{(j)} \sim
  \mathrm{Bernoulli}(\hat{\omega}),
  \qquad
  i=1,\ldots,n, \quad j=1,\ldots,N,
\end{equation}
%
where \(n\) and \(N\) denote the number of patients and the total number of bootstrap replicates, respectively, and $j$ is the index of the bootstrapped sample. Biomarker values were then drawn conditionally as
%
\[
  X_i^{(j)} \mid Y_i^{(j)} =\left \{
  \begin{array}{lll}
    \mathrm{BetaPrime}(\hat{a},\hat{b},\hat{s})&: & Y_i^{(j)}=1, \\
    \mathrm{Gamma}(\hat{k},\hat{\vartheta})&: & Y_i^{(j)}=0.
  \end{array} \right.
\]
%
Each simulated dataset was refitted to obtain
%
\[
  \hat{\bTheta}^{(j)}, \qquad x_{50}^{(j)}.
\]
%
Confidence intervals were again computed from empirical quantiles across bootstrap replicates. This procedure quantifies uncertainty under the fitted generative model structure.

\paragraph{Primary inferential target.}

Although we computed confidence bands for the fitted risk function, the principal inferential quantity was the midpoint \(x_{50}\). Accordingly, we obtained confidence intervals for \(x_{50}\) directly from its propagated or resampled distribution, rather than indirectly from the fitted curve bands.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliographystyle{jphysicsB}
%\bibliographystyle{agsm}
\bibliographystyle{jphysicsB}
\bibliography{Ref}

\end{document}