Active Learning-Based Multistage Sequential Decision-Making Model with Application on Common Bile Duct Stone Evaluation

Multistage sequential decision-making is commonly seen in many real-world applications for the sake of accuracy, efficiency, and cost-effectiveness, such as heathcare audit [1], optimal investment [2], energy [3], optimization [6], treatment [5], and design of experiments [4]. Meanwhile, active learning is a special case of machine learning in which a learning algorithm can interactively query a user or some other information source to label new data points with the desired outputs [7], and has been integrated with decision-making in recent years in various contexts such as imbalanced data [8, 9] and Bayesian decision-making [10]. In this paper, we aim at developing an active learning-based multistage sequential decision-making model to assist doctors making reliable diagnostic decisions and treatment recommendations in a cost-effective and convenient manner.

Our research is motivated from the evaluation of gallstone. Gallstone are solid particles that can form from cholesterol, bilirubin, and other substances within the gallbladder. These densities are often benign when localized to the gallbladder, but can cause pain, infection, and liver damage when they become stuck in the common bile duct (CBD) and impede the flow of bile into the digestive tract. A stone that becomes impacted in the CBD can be difficult to detect definitively, but its presence requires a procedural intervention, endoscopic retrograde cholangiopancreatography (ERCP). As with most procedures, ERCP requires anesthesia and carries the risk of complications, including pancreatitis, infection, and bleeding. As a result, it is imperative to ensure that ERCP is only performed when there is a stone definitively obstructing the duct. This is particularly true in children [11]. Thus, one important objective is how to efficiently predict the presence of a CBD stone in a child with high accuracy and specificity.

From the machine learning or statistical point of view, this might be straightforward: one would like to collect as many informative features as possible for each patient, which will lead to higher prediction accuracy. However, from the medical point of view, this is highly non-trivial, since some informative features could be expensive and hard to obtain in many real-world scenarios. Obtaining additional clinical data often comes with a price, including delayed time to intervention, cost, and complications. This price can be particularly high in pediatric care. For instance, while magnetic resonance imaging (MRI) is the most accurate way to diagnose a CBD stone, only 5% of pediatric patients in our real dataset have such diagnostic, as pediatric MRI often involves sedation, long scheduling time, etc. Indeed a less cost-effective approach such as laboratory tests or ultrasound (US) often might be sufficient. Thus a common clinical challenge is how to reduce the cost of data collection for diagnosis decision-making, yet avoid computational issues associated with limited or incomplete data.

In this article, a multistage sequential decision-making model is developed to assist doctors making reliable diagnostic decisions and treatment recommendations with high accuracy in a cost-effective and convenient manner. In our approach, some patients are allowed to take fewer examinations if their symptoms are apparent according to early-stage features, and some patients will be recommended to take more examinations if their symptoms are ambiguous and doctors have difficulty making reliable treatment recommendations.

There are several benefits to the proposed model. This novel platform allows the healthcare provider to actively collect only the necessary data sequentially rather than collecting all the data at inception. As a result, the diagnostic process can be faster, more convenient, and cost effective, while retaining significant accuracy. In the research setting, our model offers a novel solution to address challenges raised by incomplete data through offering sequential decision-making in an active learning way. The active learning-based model allows some missing data yet guarantees reliability and prediction accuracy, which makes the incomplete data missing “at control” rather than missing at random.

In the proposed framework, one important thing is to identify a proper predictive model for each stage as the building block. General categorization is required to utilize the model. Patients must be stratified to: certainly healthy (“0”), indeterminate (“0.5”), and certainly sick (“1”). After comparing the performance of different machine learning models, including neural networks, random forest, logistic regression, we selected the logistic regression model, which provides stability and accuracy with easy interpretability, see our previous work [11]. In the sections to follow, a generalized logistic regression model will be utilized as a building block for each stage.

Multiclass logistic regression, also called softmax/multinoulli/multinomial regression, is a classification method that generalizes logistic regression to multiclass problems ($K>2$). The idea behind this classification model is one-vs-all or one-vs-rest, which means $K$ binary classifier models will be trained to ensure every class has a binary classifier to identify whether a sample belongs to this specific class or not. This model enables multiple outcomes but does not involve the ordering of the categorical dependent variable. The ordinal logistic regression model, also called the ordered logit model or proportional odds model, overcomes this limitation by calculating the log of the odds of cumulative events [14, 15, 16]. Unlike simple logistic regression, ordinal logistic models consider the probability of an event and all the events that are below the focal event in the ordered hierarchy. The existing ordinal logistic regression model satisfies our needs in single stages but failed to deal with multistage scenarios simultaneously. Therefore, in this research, we’ll develop our own multistage model with the existing ordinal logistic regression model as a building block for each stage, but estimate the parameters in all stages together.

The contribution of our work is twofold. In terms of statistical analysis, an active learning model is proposed in a multistage decision-making process. This learning algorithm can interactively query a patient collect additional clinical data for collecting necessary information [17]. Also, the proposed model introduces a novel way of dealing with incomplete data by allowing some data to be absent and modeling with incomplete but necessary data directly, rather than guessing the incomplete data with uncertainty. In terms of healthcare application, a general strategy is developed to help the doctor making diagnostic decisions and treatment recommendations sequentially when some laboratory data are more expensive, invasive, or time-consuming to collect than other data. By recommending requisite testing to those patients in whom symptoms do not support a definitive diagnosis, the proposed process is fast, cost-effective, yet reliable. Also, the methodology applies to a large variety of disease diagnose that involve a sequence of indirect examinations and some unpleasant or inaccurate examinations.

Note that our proposed model/method is closely related to incomplete data analysis in the statistics literature, where many widely-used methods have been developed, e.g., simple deletion on feature or record, mean/median/mode/zeros substitution, regression imputation, last observation carried forward, EM, etc. [13, 12]. However, these existing methods were developed from two primary approaches. One is deleting the records or features with missing values, which will lead to a waste of available data and some intrinsic problems associated with limited data size. The other approach is “guessing” or “learning” the missing data from the observed data, which can be dangerous for clinical modeling since the synthetic data brings extra uncertainties. Our approach is different from these existing methods in the sense that we actively decide which kinds of data to be observed and which kinds of data to be missing.

In Section 2, the proposed model is described. Section 3 discusses the novel approach to parameter estimation for the proposed model. In Section 4, the performance of the proposed methodology on synthetic data is reported. Finally, in section 5, the analytical result on real clinical data is demonstrated.

In this section, we first present a high-level framework of our proposed active learning-based multistage sequential decision-making model, and then introduce necessary notation to represent observed data in our framework. Finally, our model is rigorously developed based on the existing ordinal logistic regression model.

At the high-level, we propose to use the maximum likelihood estimation (MLE) method to estimate not only the unknown $\beta$’s values in different stages, but also the cutoff values $\left\{L_{1}^{*},U_{1}^{*},C_{2}^{*}\right\}$. Note that the coefficients $\beta$’s can also be estimated by the penalized likelihood estimation framework such as LASSO for variable selection when there are a large of available features, which is beyond the scope of this article and will be presented elsewhere. Here we assume that we pre-screen all features, and only use those useful features in our model and medical decision making.

To better present our main ideas, let us begin with a single stage with a ordinal logistic regression model. To simplify the notation, we assume that the response $Y$ results from the underlying latent variable $Y^{*}$ by using cutoff points $\theta_{1}<\theta_{2}$:

where $\theta_{0}=-\infty$ and $\theta_{3}=\infty$ are introduced for the technical reason to simplify the notation. Also note that at the last stage, we have $\theta_{1}=\theta_{2}$, which will force one to make a decision $\{Y=0\}$ or $\{Y=1\}.$ Under this notation, the conditional distribution of $Y$ is given by

where $\boldsymbol{\beta},\theta_{j}$ are model coefficients and cutoff points to be estimated from the data, $P$ is the abbreviation of probability, and $F$ is the CDF of the error term $\boldsymbol{\varepsilon}$ which takes on the role of the inverse link function.

For a given stage or a given ordinal logistic regression model with $N$ subjects, the joint log-likelihood function becomes

It has been well investigated to maximize this log-likelihood objective function which possesses several nice properties. McCullagh presented a Fisher scoring algorithm for ML estimation, expressing the likelihood function using cumulative probabilities. McCullagh showed that a sufficiently large $n$ guarantees a unique maximum of the likelihood [14, 19]. Burridge and Pratt showed that the log-likelihood is concave for many cumulative link models, including the logit, probit, and complementary log-log [19]. Iterative algorithms usually converge rapidly to the ML estimates. The optimization computing algorithms are developed based on these useful properties. The existing approaches compute the MLE by iteratively reweighted least squares (IRLS). Actually, for packages in $R$ software, all GLM are fit using the IRLS estimation method for the log-likelihood.

Now let us go back to the parameter estimation of our proposed multistage model. For simplicity, denote $[\boldsymbol{\beta},\boldsymbol{\theta}]$ as $\boldsymbol{\zeta}$, and $\log\mathcal{L}\left(\boldsymbol{\zeta}\mid\boldsymbol{X},Y\right)$ as $\ell(\boldsymbol{\zeta})$. Under our notation, the joint log-likelihood function for all samples in all stages will be as follows:

where $K$ is the number of stages in the proposed model, $N_{k}$ is the number of patients in the $k^{th}$ stage, and $J_{k}$ is the number of categories in the $k^{th}$ stage for which we have

Given the joint likelihood function $\ell(\boldsymbol{\zeta})$, the next step is to simultaneously estimate $\boldsymbol{\zeta}$, all the parameters in the $K$ stages, that maximizes the likelihood function. Here we propose to adopt the gradient descent algorithm and a key step is to compute the gradient for each parameter, which can be written as:

where $f(t)=F^{\prime}(t)$ is the probability density function (PDF) of $\boldsymbol{\varepsilon}$: $f(t)=\frac{1}{(e^{t/2}+e^{-t/2})^{2}}$.

We should also remark some technical details in computation. Note that the coefficient in stage $k$, $\boldsymbol{\beta}^{(k)}$, contains $\boldsymbol{\beta}^{(k-1)}$ and newly added coefficients for features in the $k^{th}$ stage. This makes the coefficients in different stages of different lengths. To address this issue, we propose finding the gradient in each stage $\frac{\partial\ell(\boldsymbol{\zeta})}{\partial\boldsymbol{\beta}^{(k)}}$ individually first. Then, adding them together with the front ends aligned and pad the shorter ends with 0’s: $\frac{\partial\ell(\boldsymbol{\zeta})}{\partial\boldsymbol{\beta}}=\left[\frac{\partial\ell(\boldsymbol{\zeta})}{\partial\boldsymbol{\beta}^{(1)}};0,0,...,0\right]+...+\left[\frac{\partial\ell(\boldsymbol{\zeta})}{\partial\boldsymbol{\beta}^{(K)}}\right]$. In addition, the indicator function $I(y_{i}^{(j)}=j)$ can be replaced with a mask matrix with boolean values. Alternatively, the indicator function can be avoided by data preprocessing. For instance, we can split the patients into different groups to mark the identity, e.g., using $G_{k,j}$ to denote the group of patients who have label $j$ in stage $k$, see Appendix A for more details.

To numerically solve our proposed joint-optimization problem, we follow the framework of the existing $polr()$ function in the package MASS in R software, and adopt the $optim()$ function as optimization solver. The $optim()$ is a function designed for general-purpose optimization based on Nelder–Mead, quasi-Newton and conjugate-gradient algorithms. In $optim()$, the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) method, an iterative method for solving unconstrained nonlinear optimization problems [20], is adopted, since it works best with analytic gradients based on our extensive experiences.

In this section, a simulation is conducted to verify the effectiveness and efficiency of the proposed methodology. First, a 2-stage synthetic dataset, $\{X^{(1)},Y^{(1)};X^{(2)},Y^{(2)}\},$ is generated with prescribed coefficients $\boldsymbol{\beta}^{(1)},\boldsymbol{\beta}^{(2)}$ and classification thresholds $\{L^{*}_{1},U^{*}_{1},C^{*}_{2}\}$. Then, we conduct both a baseline approach and the proposed approach on the synthetic dataset to demonstrate the performance of the proposed method in the aspect of prediction accuracy, stability, interpretability, and relative efficiency.

Let $N_{1}=10,000$ denote the total number of synthetic patients at the outset. We design the features of patients as follows: $X_{1}\sim Bernoulli(0.3)$, $X_{2}\sim N(-1,1)$, $X_{3}\sim N(1,1)$, $X_{4}\sim N(0,2)$, $X_{5}\sim Bernoulli(0.4)$, $X_{6}\sim N(-1,1)$, $X_{7}\sim N(0,1)$. For stage 1, we design feature set as: $\boldsymbol{X}^{(1)}=\left\{X_{1},X_{2},X_{3},X_{4}\right\}$. For stage 2, we add 3 new features: $\boldsymbol{X}^{(2)}=\left\{X_{1},X_{2},X_{3},X_{4};X_{5},X_{6},X_{7}\right\}=\left\{\boldsymbol{X}^{(1)};X_{6},X_{7},X_{8}\right\}$. The error term $\boldsymbol{\varepsilon}$ follows standard logistic regression for both stages: $\boldsymbol{\varepsilon^{(1)}}\sim Logistic(0,1),\quad\boldsymbol{\varepsilon^{(2)}}\sim Logistic(0,1)$. Note that if the error term $\boldsymbol{\varepsilon}$ follows logistic regression with other parameters, the estimation of $\boldsymbol{\beta}$ and $\boldsymbol{\theta}$ will be scaled, see Eq. (2). The coefficients are designed as: $\boldsymbol{\beta}^{(1)}=[2,2,2,2],\quad\boldsymbol{\beta}^{(2)}=[\boldsymbol{\beta}^{(1)},4,4,4]=[2,2,2,2,4,4,4]$ which satisfys the common coefficients assumption. Next, we can calculate the latent variable $Y^{*}$ according to Eq. (1). To cut the $Y^{*}$ into ordered categories, we design the cutoff points as $\boldsymbol{\theta}=\{L_{1}^{*},U_{1}^{*},C_{2}^{*}\}=\{-2.2,2.2,0.5\}$. The value $-2.2$ and $2.2$ are chosen since $U_{1}=F(2.2)=0.9$ and $L_{1}=F(-2.2)=0.1$, where $0.9$ and $0.1$ are the designed cutoff points on probability $Y$. Via such design, we have $N_{2}\approx\frac{1}{3}N_{1}$. Given $\boldsymbol{\theta}$, the corresponding true label $Y$ can be obtained according to Eq. (2).

In the baseline approach, the 2 stages proceed independently: in stage 1, a standard ordinal logistic regression is applied on $\{\boldsymbol{Y}^{(1)},\boldsymbol{X}^{(1)}\}$ via $polr()$ function in R; and in stage 2, a standard logistic regression is applied on $\{\boldsymbol{Y}^{(2)},\boldsymbol{X}^{(2)}\}$ via the $glm()$ function in R. While in the proposed method, all the parameters in 2 stages are estimated together via MLE. To provide reliable statistical inference, we adopt the Monte-Carlo method and run the data generation and parameter estimation 100 times.

The comparison of the mean coefficients estimated from the two methods is visualized in Fig.  2. Clearly, the mean coefficients estimated from the two methods are similar and both coincide well with ground truth, indicating that both methods have little or no bias on mean parameter estimation, proving the effectiveness of both methods on estimating. To quantify the subtle difference, we can utilize mean square error (MSE) to measure the distance between the estimated coefficients and the ground truth. For coefficients in stage 1, the MSE of the estimation from the baseline method and the proposed method are $7.62\times 10^{-6}$ and $2.25\times 10^{-6}$, respectively. Meanwhile, for stage 2, the MSE of the baseline method and the proposed method are 0.0011 and 0.0003, respectively. This indicates a slight advantage of the proposed methods on mean coefficients estimation. See Appendix B for more comparison on mean prediction metrics.

Fig. 4 shows the standard deviation (std) of coefficients estimation from the two methods. It can be seen that the standard deviation of estimation from the proposed method is significantly smaller than those from baseline method in both stages. To quantify this advantage, the relative efficiency of the proposed method with respect to the baseline method, defined as $e(T_{1},T_{2})=\frac{var(T_{2})}{var(T_{1})}$, is plotted in Fig.  4, where $T_{1}$ and $T_{2}$ represent the estimators for $\hat{\boldsymbol{\beta}}$ in the proposed method and baseline method, respectively. It can be seen that the relative efficiency of the estimator in the proposed method for the coefficients in stage 1 is very big, as high as 19.38. And for the estimator for the coefficients in stage 2, the lowest relative efficiency is still $1.62$. This indicates that the efficiency of the parameter estimation is improved $62\%~{}1838\%$ in the proposed methodology. The reason for such improvement is straightforward: in the proposed methods, all available data is involved in the parameter estimation, while in the baseline method, only a small subset of the data is involved in later stages.

For the estimation of cutoff points, the baseline method and the proposed method give the mean estimations as $\{-2.201,2.192,0.536\}$ and $\{-2.199,2.197,0.506\}$, respectively. And the standard of deviation of the estimations by the baseline and proposed methods are $\{0.061,0.064,0.178\}$ and $\{0.057,0.062,0.119\}$, respectively. Compared with the ground truth, $\{\boldsymbol{\theta}\}_{ground~{}truth}=\{-2.2,2.2,0.5\}$, both methods are functional, yet the proposed approach yields a lower mean deviation and standard deviation.

Overall, the proposed 2-stage model provides a more accurate and efficient estimation since it involves all data in the parameter estimation. Also, in implementation, the proposed algorithm enjoys an advantage with regard to computational convenience since we estimated all parameters in $\zeta$ at once. Additionally, the common coefficients assumption in the proposed method reduces the number of parameters to be estimated from $4+7+2+1=14$ to $7+2+1=10$ since the coefficients for features in stage 1 are not repeatedly estimated in later stages, but instead assumed to be consistent.

In this section, we analyze a real data set in [11] and demonstrate the usefulness of our proposed approach in the sense of being able to learn from data to effectively evaluate gallstone. In particular, the real clinical data is modeled as a 2-stage sequential decision-making model. We have $N_{1}=316$ patients in the first stage. Based on the missing profile, we filter the patients who have features in category 2 into the second stage and label them as $Y_{i}^{(1)}=0.5,i=1,...,N_{2}$. As a result, $N_{2}=189$. Note that here we made the assumption that if a patient has available data in the next stage, the patient should be labeled as “0.5” at the current stage. The reason for such assumption is that we believe that the provider only recommends further diagnostic interventions to a patient who has indeterminate presentation at the current stage (“0.5”). In other words, the real data comes from the provider’s knowledge and experience, and the missing pattern is not missing at random but deliberately interfered with. Note that features in $\boldsymbol{X}^{(1)}$ include hemolytic disease, gender, age, body mass index (BMI), total bilirubin, alanine transaminase (ALT), aspartate transaminase (AST), lipase, amalyse, gamma-glutamyl transferase (GGT), alkaline phosphatase. These features include demographic and laboratory tests. For features in $\boldsymbol{X}^{(2)}$, based on $\boldsymbol{X}^{(1)}$, it involves imaging data, including CBD diameter and the presence of CBDS by ultrasound (US) and magnetic resonance imaging (MRI) examinations, respectively. See Table 1 in Appendix C for descriptive statistics.

Fig.  5 visualizes the coefficients for features in category 1 estimated by baseline and proposed approach, respectively. It can be seen that with the baseline approach, the coefficients estimated by the two individual models are similar for most features, $X_{2}\sim X_{11}$, but significantly different for $X_{1}$. In the proposed 2-stage model, the coefficients estimated for feature $X_{1}\sim X_{11}$ seem to resemble a weighted average of the coefficients obtained from baseline individual models. Also, compared with the baseline method, the coefficients estimated by the proposed method have no dramatic spikes. Also, with the common coefficients assumption, newly added features will not affect the weights of existing ones. More importantly, in the proposed method, no matter how many stages a provider utilized, the calculation can be done at once efficiently.

Our model and method can be useful for providing the treatment recommendation for a new patient, and a diagnostic mobile App was developed in 2021 in Apple Store based on this research. To be concrete, when a new patient presents, given the trained proposed method, the provider (e.g., clinician or doctor), can first recommend the patient undergo diagnostic tests included in category 1, and calculate the $\hat{Y}^{(1)}$ accordingly. Based on the $\hat{L}_{1}$ and $\hat{U}_{1}$, the patient will be labeled as healthy (“1”), sick (“0”), or indeterminate and proceed to the next stage (“0.5”). Then, only when the patient’s presentation is indeterminate, additional testing such as ultrasound will be undertaken. In other words, only the necessary data are actively collected sequentially for sequential decision-making. This approach simplifies and expedites the diagnostic workup. Additionally, the provider can adjust the multistage sequential decision-making model flexibly by designing the number of stages and the features in each stage, also see Section 2 for the details on the feature grouping recommendations.

In this study, to augment healthcare clinical decision-making based on limited examinations, a multistage sequential decision-making model is developed to actively collect only minimal necessary diagnostic data.To illustrate the proposed method, we optimized a 2-stage model on synthetic data in the simulation study and applied the proposed algorithm to pediatric patients with suspicion for obstructive common bile duct stones in an effort to more efficiently and accurately arrive at a diagnosis that allowed expedited intervention with ERCP. In both the simulation study and real case study on collected pediatric data, the proposed method is more effective, stable, interpretable, and simpler in estimating the coefficients and cutoff points for all stages. It estimates all parameters simultaneously and fully considers the cumulative characteristics between these stages by making common coefficients assumption.

One limitation of the proposed methodology is that it involves all variables without selection. A potential future improvement to the model will be the addition of $L_{1}$ regularization term to the objective function to realize the variable selection.

Also, we can extend our proposed algorithm to satisfy more real-world cases with special requirements. For example, consideration of expense, delay to intervention, and need for invasive testing can be specified and minimized with modification to the objective function as $\min_{N_{1},N_{2},\zeta}-\log\mathcal{L}(\zeta\mid X,Y)+\lambda\left(c_{1}N_{1}+c_{2}N_{2}\right)$, where the $c_{1},c_{2}$ are the relevant cost per patient in stage 1 and stage 2, respectively. $\lambda$ is the tuning parameter to adjust the weight between the negative likelihood and the relevant cost to the test. Note that here the $N_{1},N_{2}$ can also be variables to be optimized over since the relevant cost depends on the number of patients in each stage.

In healthcare settings, the penalty for misdiagnosis can vary based on the patient and diagnosis. For example, false-positive detection may lead to unnecessary procedures or excess intervention, while false-negative detection may result in missed or delayed treatment. For such cases with weighted penalties on different types of misclassification, we can modify the objective function by specifying the penalties as $\min_{\zeta}-\log\mathcal{L}(\zeta\mid X,Y)-\lambda_{1}\sum_{i}I\left(\left(y_{i}=0.5\mid 1\right),\hat{y}_{i}=0\right)-\lambda_{2}\sum_{i}I\left(\left(y_{i}=0\mid 0.5\right)\&\hat{y}_{i}=1\right)$, where the $\lambda_{1}$ and $\lambda_{2}$ are the tuning parameter for false-negative prediction and false-positive prediction, respectively.

The proposed algorithm is also applicable to a broad range of decision-making scenarios outside of healthcare as it utilizes active learning to drive the sequential collection of minimal requisite information, rather than all possible information simultaneously. This will be an area of future development.