Simultaneous Modeling of Disease Screening and Severity Prediction: A Multi-task and Sparse Regularization Approach

One of the important issues in biomedical research is the exploration of biomarkers and the construction of clinical prediction models. Biomarkers are physiological indicators or biomolecules used for various purposes such as disease screening, prognosis, and estimation of treatment effects, and they are widely used in clinical practice (Biomarkers Definitions Working Group., 2001). This paper focuses specifically on biomarkers used for disease screening, aiming to propose models and estimation methods for the exploration of biomarkers that are not only useful for screening but can also predict the disease severity at screening. First, we will show examples of such muti-purpose biomarkers.

Alpha-fetoprotein (AFP), Des-gamma-carboxyprothrombin, and Lens culinaris agglutinin-reactive fraction of AFP are well-known biomarkers for hepatocellular carcinoma. These biomarkers are measured through blood testing and serve as less invasive diagnostic markers for HCC. Furthermore, it has been reported that these markers are associated with prognosis and cancer stage (Omata et al., 2017; for the Study of the Liver, 2018). Another example is neuroblastoma, which is a type of childhood cancer known to have a wide range of severity (Irwin et al., 2021). It is known that neuroblastoma patients excrete specific metabolites in their urine, and conventionally, molecules such as homovanillic acid (HVA) and vanillylmandelic acid (VMA) have been used for disease screening. However, while these are highly effective in screening neuroblastoma, their association with severity is weak. Recently, a metabolomics study has revealed some markers are useful for screening and related to severity, beyond HVA and VMA (Amano et al., 2024).

How should we take into account disease screening and severity prediction at one time? A simple way would be to stack them into one ordinal variable. In other words, to define the absence of disease as the lowest category and set higher categories for increasing disease severity, treating it as an ordinal categorical variable. In regression analysis, an ordinal categorical variable is one of the most common types of response variables. A well-known regression model for ordinal responses is the cumulative logit model (CLM, a.k.a. the proportional odds model; McCullagh, 1980). CLM assumes a continuous latent variable behind the response variable, and the response variable level increases when the latent variable exceeds certain thresholds. CLM is also widely used in biomedical studies (e.g., Muchie, 2016; Jani et al., 2017). However, this approach may have an issue applying to our problem; CLM assumes an identical relationship between the response and the predictor across all levels of the response variable, which is referred to as parallelism assumption. The parallelism assumption is questionable for the combined response because the two tasks, screening and severity prediction, do not necessarily share the same structure; in other words, the same predictors do not necessarily contribute in the same way to screening and severity prediction. One possible solution to this issue is to relax the parallelism assumption using a varying coefficient version of CLM, which we call non-parallel CLM (NPCLM). However, as discussed in Section 2, the varying coefficient approach could make the estimation and interpretation of the model difficult.

To address the problem of modeling the combined task, we propose a model incorporating the idea of multi-task learning, we aim to effectively construct a prediction model and discover beneficial biomarkers, even in high-dimensional settings. Multi-task learning is a machine learning approach where multiple models are trained simultaneously, leveraging shared information and structures across the tasks to improve predictive performance (Caruana, 1997). Multi-task learning is applied in many fields such as computer vision, natural language processing, and web applications (Zhang and Yang, 2017). Furthermore, multi-task learning is being incorporated in the area of biomarker discovery. For example, Wang et al. (2020) introduced multitask learning to identify genes expressed across multiple types of cancer. While there are not many examples of applying multi-task learning to ordinal regression, Xiao et al. (2023), for instance, deals with multi-task learning for parallel ordinal regression problems.

This paper is organized as follows. In Section 2, we review CLM, NPCLM, and related extensions. In Section 3, we propose a novel prediction model and also discuss its relationships with other categorical models and potential extensions. In Section 4, we introduce an estimation method using sparse regularization and an optimization algorithm. The proposed models employ structural sparse penalties to exploit the common structure between screening and severity prediction. In Section 5, we present the results of simulation experiments, considering various structures between screening and severity prediction, and clarifying under what conditions the proposed method performs well. In Section 6, we report the results of real data analysis. Finally, we make some concluding discussion in Section 7.

For ordinal responses, CLM is one of the most popular regression models (McCullagh, 1980; Agresti, 2010). Let

$\displaystyle Y_{i}\in\{0,1,\ldots,K\}~{}(:=\mathcal{Y}),~{}~{}X_{i}\in\mathcal{X},~{}~{}~{}i=1,\ldots,n,$

the response and the predictors, where $\mathcal{X}$ is the $p$-dimensional predictor space, and the set $\mathcal{Y}$ is equipped with an ordinal relation of natural numbers. CLM is defined as

$\displaystyle\mathrm{logit}~{}\mathbb{P}(Y_{i}\leq k\mid X_{i})=\alpha_{k}+X_{i}^{T}{\boldsymbol{\beta}},~{}~{}~{}k=1,\ldots,K-1,$

where $\alpha_{k}\in\mathbb{R}$ and $\boldsymbol{\beta}\in\mathbb{R}^{p}$ are the intercepts and regression coefficients, and $\mathrm{logit}(u)=\log\{u/(1-u)\}~{}~{}(u\in(0,1))$. CLM has an important interpretation. Suppose that there is a latent continuous variable $Y^{*}\in\mathbb{R}$ behind the response and that $Y^{*}_{i}$ is associated with the predictors as

$$\displaystyle Y^{*}_{i}=X_{i}^{T}{\tilde{\boldsymbol{\beta}}}+\varepsilon_{i},~{}~{}~{}~{}\varepsilon_{i}\sim Logistic(0,1),$$

where $\tilde{\boldsymbol{\beta}}\in\mathbb{R}^{p}$ is the regression parameters for $Y^{*}_{i}$, and $Logistic(0,1)$ is the standard logistic distribution. Then, by defining $Y_{i}=k~{}\text{iff}~{}\tilde{\alpha}_{k-1}<Y^{*}_{i}\leq\tilde{\alpha}_{k}$ with an increasing sequence $\{\tilde{\alpha}_{k}\}_{k=0}^{K}$ with $\tilde{\alpha}_{0}=-\infty$ and $\tilde{\alpha}_{K}=\infty$, we have

$$\displaystyle\mathrm{logit}~{}\mathbb{P}(Y_{i}\leq k\mid X_{i})=\mathrm{logit}~{}\mathbb{P}(Y_{i}^{*}\leq\tilde{\alpha}_{k}\mid X_{i})=\mathrm{logit}~{}F_{\varepsilon}(\tilde{\alpha}_{k}-X_{i}^{T}{\tilde{\boldsymbol{\beta}}})=\tilde{\alpha}_{k}-X_{i}^{T}{\tilde{\boldsymbol{\beta}}},$$

where $F_{\varepsilon}$ is the cumulative distribution function of $Logistic(0,1)$. By replacing $\alpha_{k}=\tilde{\alpha}_{k}$ for all $k\in\{1,\ldots,K-1\}$ and $\boldsymbol{\beta}=-\tilde{\boldsymbol{\beta}}$, we can see that this model is equivalent to CLM. Therefore, we can see $\alpha_{k}$ as the thresholds determining the class $Y_{i}$ based on the latent variable $Y_{i}^{*}$. The linear functions $\alpha_{k}+X_{i}^{T}\boldsymbol{\beta}$ representing the log-odds of the cumulative probability of each level are all parallel because they share the same slope $\boldsymbol{\beta}$. The parallelism assumption is necessary to ensure that the model is valid in the sense that the conditional cumulative probability $\mathbb{P}(Y_{i}\leq k\mid X_{i}=x)$ derived from CLM is monotonic at any given $x\in\mathcal{X}$ (e.g., Okuno and Harada (2024)).

The non-parallel CLM (NPCLM), sometimes called the non-proportional odds model, is the model with different slopes for each response level. This implies the conditional cumulative probability curves of NPCLM can be non-monotone. Such curves violate the appropriate ordering of the response probabilities. So NPCLM is more flexible and expressive than CLM, but it is valid on some restricted subspace of $\mathcal{X}$. Peterson et al. (1990) proposed a model that is intermediate between CLM and NPCLM (Peterson and Harrell, 1990). This partial proportional odds model is defined as follows:

$$\displaystyle\mathrm{logit}~{}\mathbb{P}(Y_{i}\leq k\mid X_{i})=\alpha_{k}+X_{i}^{T}{\boldsymbol{\beta}}+X_{i}^{T}{\boldsymbol{\gamma}}_{k},~{}~{}~{}k=1,\ldots,K-1,$$

(1)

where $\boldsymbol{\gamma}_{k}$ is the variable slope of the $k$th level. This model is designed to capture the homogeneous effect by $\boldsymbol{\beta}$ and the heterogeneous effect by $\boldsymbol{\gamma}_{k}$. Still, similarly to NPCLM, this model is generally valid as a probability model only on the restricted subspace of $\mathcal{X}$. To control the balance of flexibility and monotonicity, various estimation techniques with regularization have been proposed. Wurm et al. (2021) use L1 and/or L2 norm on coefficient parameters $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}_{k}$ of the partial proportional odds model (Wurm et al., 2021). If the parallel assumption holds, then $\boldsymbol{\gamma}_{k}$ should be zero due to L1 penalization, and even when the parallel assumption is violated, the variability of $\boldsymbol{\gamma}_{k}$ is controlled by the penalties. Wurm et al. (2021) have also proposed an efficient coordinate descent algorithm (Wurm et al., 2021), which is similar to that of the lasso regression (Tibshirani, 1996; Hastie et al., 2015). Tutz et al. (2016) have introduced the penalization on the difference of the adjacent regression coefficient (Tutz and Gertheiss, 2016). For NPCLM, their penalty term is $\sum_{k}\|\boldsymbol{\gamma}_{k+1}-\boldsymbol{\gamma}_{k}\|_{2}^{2}$, resulting in smoothed coefficients between the adjacent response levels. They use the L1 penalty for the adjacent coefficients instead of the L2 penalty. With the L1 penalty, unlike the L2 penalty, the adjacent regression coefficients are expected to be estimated as the same if the parallelism assumption holds. Additionally, they have proposed an algorithm based on the Alternating Direction Method of Multipliers (ADMM; Boyd et al. (2011)).

In our setting, the non-parallel models may be helpful, but as previously discussed, they can be non-monotone for given $x$. Furthermore, although these models are linear, they are sometimes not easy to interpret. If the estimated coefficients indicate that $\mathrm{logit}~{}\mathbb{P}(Y_{i}=healthy\mid X_{i}=x)$ and $\mathrm{logit}~{}\mathbb{P}(Y_{i}\leq mild\mid X_{i}=x)$ have different coefficients with opposite sign for $X$, it may be difficult to understand, given that the event $\{Y_{i}\leq mild,X_{i}=x\}$ includes $\{Y_{i}=healthy,X_{i}=x\}$. In the next section, we propose a multi-task learning approach that maintains monotonicity over the entire $\mathcal{X}$ and possesses a flexibility that allows it to capture the different structures between the tasks of screening and severity prediction.

As noted in Section 3, we cannot exploit the shared structure of the two tasks, screening and severity prediction, through simple maximum likelihood estimation. Therefore, we employ a penalized maximum likelihood approach.

In this section, we perform numerical experiments to compare the proposed and existing methods in prediction and variable selection performance in several scenarios. In all scenarios, the sample size is set to 300. Response $Y$ has four levels $\{0,1,2,3\}$, generated by thresholding the latent variable $Y^{*}$ and $Y^{**}$. Predictors are generated from standard normal distributions that are independent of each other. The dimensions of the predictors are set to four values for each scenario: 75, 150, 300, and 450. As described in Table 1, five scenarios with different true structures, in which 10 to 20 predictors are truly relevant to the response, are set up to account for various situations.

Details on the simulation models are found in Appendix C. Note that Scenario 2 appears to satisfy the parallelism assumption, but not as mentioned in Section 3.2. Scenarios 2 to 5 do not meet the parallelism assumption, and in particular, Scenarios 3 and 5 are in serious violation.

The proposed and existing methods include the MtCLM with three types of penalties (L1, L1 + Group lasso, and L1 + Fused lasso), the L1-penalized logistic regression (glmnet package; only for screening), and the L1-penalized parallel and non-parallel CLM (ordinalNet package (Wurm et al., 2021)). All these methods have tuning parameters. For the proposed methods, $\lambda_{11}$ is chosen from $(0.01,0.05)$, $\lambda_{12}$ is chosen from $(0.05,0.1)$, and $\lambda_{F},\lambda_{G}$ are chosen from $(0,0.01,0.05)$ using 5-fold CV. Similarly, the tuning parameters of the existing methods are chosen from $(0.01,0.05,0.1)$ based on a 5-fold CV. The screening performance is evaluated using the Area Under the Curve of Receiver Operating Characteristic (ROC-AUC) based on the prediction probabilities and the F1 score when classified with a threshold of 0.5 for the prediction probabilities. The performance in severity prediction, including screening, is considered as the performance in the classification of four categories (0, 1, 2, 3), and is evaluated using Accuracy, Mean Absolute Error (MAE), and Kendall’s $\tau$, which is a measure of concordance for ordinal categories. The performance in variable selection is evaluated in terms of the power and false discovery rate (FDR) based on the ground truth. These metrics are calculated for the three types of MtCLM based on the number of regression coefficients estimated as non-zero and the number of truly non-zero regression coefficients, for scenarios where the true structure has a hierarchical structure of MtCLM. These measures are calculated over 30 simulations for each scenario. Details about the simulation settings and the evaluation measures are found in Appendix C.

Figure 2 shows the screening performance under the situation with high-dimensional predictors ($N=300,p=450$). The penalized logistic regression, which did not model the severity, showed similar performance across all scenarios and metrics. The penalized CLM exhibited superior performance, particularly in AUC, in Scenario 1, where the parallelism assumption holds, and Scenario 2, where the situation is close to Scenario 1. This result implies that, under appropriate conditions, considering severity can improve screening performance even when tackling a screening task. On the other hand, the NPCLM performed worse than the other methods in all scenarios. This performance degradation can also be attributed to the defect in non-parallel models described in Section 2. Comparing the three MtCLMs, there was no significant difference in screening performance regardless of the penalty used, but the MtCLM (L1) showed slightly larger variability in performance compared to the others. The MtCLM consistently demonstrated performance comparable to or slightly lower than the penalized logistic regression across all scenarios. By making the model more flexible compared to CLM, the performance remained stable even in scenarios where the two tasks had almost no similar structure. Furthermore, in the scenarios with high commonality between tasks (Scenarios 1, 2, and 4), it showed performance comparable to the penalized logistic regression, suggesting that it effectively leveraged the shared structure in appropriate situations despite increased variance due to increased parameters. In the settings with lower dimensional predictors, we obtained similar results; however, it is worth mentioning that the MtCLM slightly outperformed the penalized logistic regression in some scenarios when $p=75$ (See Figure A2 in Appendix Appendix D).

Predictive Performance for the the ordinal responses, namely, the joint task of screening and severity prediction is presented in Figure 3.

As with screening, in Scenarios 1 and 2, where the structure of both tasks is highly similar, the penalized CLM demonstrated high performance. Particularly in terms of Accuracy, the CLM exceeded the others; however, in terms of concordance and MAE, MtCLM (L1 + Fused) was comparable. The penalized NPCLM failed to perform well in all scenarios, and the MtCLM consistently showed high performance across all scenarios, especially in Scenarios 1 and 2 the Fused Lasso penalty worked effectively. The large variance of MtCLM (L1) was also similar to the screening performance results. Besides, MtCLM (L1 + Group) slightly outperformed others in terms of Accuracy in scenarios 3, 4, and 5. When the predictors are of lower dimensions ($p=75$), MtCLM (particularly L1 + Fused) demonstrated performance equal to or better than CLM in Scenarios 1 and 2. Additionally, in cases where the same variables are related to the response, though in different directions, as in Scenario 3, MtCLM (L1 + Group) performed relatively well.

In variable selection, the results were compared among the variants of MtCLM (Figure 4). Overall, as the dimensionality of the predictors increased, the number of false discoveries also increased. In our settings, all methods were able to select most variables that are truly associated with the response. The FDR varied depending on the choice of penalty term. When using the fused lasso penalty, especially in scenarios 3, 4, and 5, the FDR was very high, and many variables unrelated to the response were selected. When using the L1 penalty alone, the FDR was lower in many scenarios on average, but the variability of the FDR was the highest. When using the group lasso penalty, in scenarios 2 and 3, it showed detection power similar to that of using the L1 penalty alone, and the FDR was relatively lower and more stable than with the fused lasso penalty.

In summary, when there is a hierarchical structure in ordinal responses, if the structures between hierarchies are highly similar, or if the parallelism assumption is believed to hold, CLM is highly effective. Moreover, even when the primary goal is the classification of the first hierarchy, namely screening, modeling the response as an ordinal response including severity may improve screning performance. However, considering the possibility of structural differences between hierarchies, a model that accounts for these differences, such as MtCLM, should be used. As seen across the five scenarios, MtCLM effectively works even in cases where the structures between hierarchies largely differ, even with penalties like the fused lasso that emphasize shared structures. Furthermore, when the structures are well-shared between hierarchies, the fused lasso penalty is particularly effective, performing comparably to CLM. The proposed methods are said to be adaptive to the data and robust to structural differences, as they improve performance by leveraging commonalities between hierarchies when structures are highly shared and fit each task as a different model when commonalities are poor.

In this study, we have discussed the combined problem of screening and severity prediction, considering that some biomarkers are associated not only with the presence of disease but also with disease severity. To address this problem, we have proposed a multi-task approach with structured sparse regularization, MtCLM, leveraging the shared structure between the two tasks. MtCLM is an ordinal regression model, which is more flexible than CLM, but unlike NPCLM, it is valid across the entire predictor space and interpretable. The performance of MtCLM on prediction and variable selection was verified by simulation, and an example of real data analysis was provided.

This study suggests the usefulness of a multi-task learning approach for tasks that are expected to have shared structures in medicine, such as screening and severity prediction of the same disease. The idea of multi-task learning is not considered to be prevalent in the field of biomedical research. However, many diseases share risk factors and genetic abnormalities (e.g., hyperlipemia is a risk factor for cardiovascular and cerebrovascular diseases (e.g., Kopin and Lowenstein (2017)), and the HER2 gene is amplified in ovarian and gastric cancers in addition to breast cancer (e.g., Gravalos and Jimeno (2008))). Therefore, the idea of multi-task learning that takes advantage of commonalities among tasks may be useful in building prediction models or investigating prognostic factors using limited data.

Kazuharu Harada is partially supported by JSPS KAKENHI Grant Number 22K21286. Shuichi Kawano is partially supported by JSPS KAKENHI Grant Numbers JP23K11008, JP23H03352, and JP23H00809. Masataka Taguri is partially supported by JSPS KAKENHI Grant Number 24K14862.

During the preparation of this work, the authors used ChatGPT 4/4o (OpenAI Inc.) in order to improve English writing. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.