The analysis of biomedical data often aims to identify a specific (small) set of characteristics or biomarkers that will allow for the most accurate and efficient discrimination between groups. For example, it is assumed that an unknown combination of characteristics offers the best possibility to distinguish a group of patients with a certain symptom or disease from all other groups. A wide variety of statistical and machine learning tools have been developed to select the optimal feature set through an analysis of labeled data (Bair and Tibshirani, 2004; Hastie et al., 2009; Luo and Ren, 2014; Taylor and Tibshirani, 2015). In practice, however, mistakes in the assignment of labels and features during the data acquisition procedure may introduce serious bias when common methods are applied or even prohibit their use. Erroneous assignments occur in many studies for a variety of reasons: experimental errors, differences in platforms used to acquire data from different groups, differences in protocols used to post-process data, or intrinsic difficulties in distinguishing the case and control groups (Lam et al., 2011; ORawe et al., 2013; Ross et al., 2013; Weißbach et al., 2021). This is particularly evident in the case of coronavirus data emerging from different sources where it has been shown that data variability is an important factor concerning the usability of such data for machine learning (Sáez et al., 2020). Furthermore, diagnoses for many diseases only reach a certain level of confidence; in some cases (e.g., Alzheimer’s disease), a 100% diagnosis is only possible during a postmortem autopsy (Gomez-Nicola and Boche, 2015). This means that subjects with a negative diagnosis might nevertheless carry a latent form of a disease without showing symptoms yet. Assigning those individuals to the control group in a (bio)medical study can introduce a strong source of errors and can certainly have severe consequences for the patients, if this mislabeling cannot be identified, and correct medical treatment will be withheld.

Most methods that are commonly applied to deal with the problem of mislabeling are based on detection of anomalies and outliers and aim to avoid the problem by thoroughly cleaning the data during preprocessing, removing those points that appear to be mislabeled (Brodley and Friedl, 1999; Barandela and Gasca, 2000; Teng, 2001; Jiang and Zhou, 2004; Frénay and Kaban, 2014; Frenay and Verleysen, 2014). They detect outliers because they deviate significantly from a model that has already been imposed on a data subset—which is again—presumed to be correctly labeled (Chandola et al., 2009). This cleaning however can be problematic because it first may not be possible if too little is known to determine what might be mislabeled, and second because it can severely reduce the size of the available data, making statistical results less reliable (Teng, 2001; Bootkrajang and Kabán, 2012). Popular methods to analyze data based on supervised (Tibshirani, 1996; Hastie et al., 2009; Luo and Ren, 2014) and semi-supervised machine learning methods (Moya and Hush, 1996; Bair and Tibshirani, 2004; Rodionova et al., 2016) for labeled data analysis (e.g., generalized linear models and neuronal networks) are also implicitly based on an assumption that at least one of the groups to be discriminated has been labeled perfectly or at least assume a subset of perfectly labeled data (Hendrycks et al., 2019). On the other side, common unsupervised methods (such as hidden Markov models, Bayesian mixture models (Frühwirth-Schnatter, 2006; Todorov et al., 2020), and advanced clustering methods (Andreopoulos et al., 2009; Gerber and Horenko, 2015; Rodrigues et al., 2021)) ignore any prior assignments of data to labels and groups in a given dataset. Herewith, however, a lot of valuable information is lost.

Also in a broader context of unsupervised data anomaly detection, two major method families like the one-class support vector machines (OCSVM, sometimes also referred to as one-class learning methods) (Choi, 2009; Zhu et al., 2016) and isolation forests (IF, anomaly detection algorithms based on random forest ideas) (Liu et al., 2008; Hariri et al., 2021) rely either on the explicit knowledge of some subsets of correctly identified data anomalies that can be used for training or on knowing the exact proportion of anomalous data in the given dataset. In the latter case, providing the exact proportion of anomalous data in OCSVM and IF allows identifying the exact value of the anomaly threshold that can be used to separate normal data from anomalous data. The primary aim of this study is providing a robust computational mislabeling inference procedure for generalized linear models (e.g., logistic regressions)—an algorithmic procedure that does neither rely on the explicit knowledge of the particular mislabeled data instances nor on the knowledge of the exact proportion of the mislabeled data in the given dataset. Instead, the introduced procedure relies on the knowledge of the upper bound for the mislabeled data proportion and deploys the tools from information theory (like Akaike information criterion) to infer the optimal logistic model and the optimal class-specific mislabeling probability matrix. We thus address the currently existing methodological gap and propose a scalable method that can realistically be applied in the analysis of biomedical data. The method permits reduced sets of discriminative features to be inferred while co-estimating the group-specific data mislabeling risks. We apply the method to two examples: breast cancer diagnoses based on radiographs (example 1) and an analysis of genomic features from a Wellderly cohort (example 2) (Erikson et al., 2016)—a cohort of people who live to be 80 years or more without having experienced a serious or chronic disease. Analysis of a synthetic dataset (i.e., a dataset created by a generalized linear model with known parameters and known group-specific mislabeling, which mimics the breast cancer imaging data from example 1) is shown in Section 3 of the Supplementary Material.

Here, we give a brief description of the methodology. Detailed mathematical derivation and an investigation of its mathematical properties (a case of mislabeled Bernoulli trials, proofs of conditions for existence and uniqueness of solutions, monotonicity, and convergence of the numerical method) can be found in Lemmas 1–5 in the Supplementary Material. We consider a problem of analyzing the labeled datasets X , Y obs that are grouped into N _g cohorts/groups, with T _g being the number of instances, for example, the number of patients in the cohort/group g.

For every data instance t , g (for every patient number t in the group g), we would like to identify a relation between a vector of features X _t,g (an n-dimensional vector containing, e.g., the genotype, the age, and some other patient-specific information) and a “true”—but directly unobserved—categorical label Y _t,g . This “true” label is taking values in the finite set of m categories y = {y ₁, y ₂, … , y _m } and represents, for example, a certain phenotype. A typical setting would be to compare features from a cohort of ill people to a control group, resulting in m = 2 representing labels like y ₁ = “sick” and y ₂ = “healthy.” We consider the “true” labels Y _t,g to be unobserved since they are not directly available. Only the observed labelings Y t , g obs are available and can be mislabeled in every instance t , g with some, yet unknown, cohort-specific mislabeling probability r i , j , g = Y t , g = y i | Y t , g obs = y j . We assume that a parametric (i.e., dependent on a vector of parameters α) discriminative model is establishing a conditional dependence between the particular “true” unobserved labeling and a particular observed feature vector X _t,g . The parametric function relating features and labels is denoted as ϕ i X t , g , α . These parametric functions can, for example, be a generalized linear model (GLM, e.g., the standard logit and probit models for m = 2) (McFadden, 1974; Friedman et al., 2010) or a neuronal network (Hastie et al., 2009).

Let y t , g , j obs = χ Y t , g obs = y j , where χ is the indicator function, taking value 1 if its argument is true and 0 otherwise. Then it can be shown that the unknown optimal model parameters α* can be inferred together with the optimal mislabeling risk matrix r* by solving the following maximization problem (see Section 1 of SI for a step-by-step derivation): α ∗ , r ∗ = argmax α , r L ( α , r ) , L ( α , r ) = ∑ g , t = 1 N g , T g 1 N g T g log ∑ i , j = 1 m y t , g , j obs r i , j , g ϕ i X t , g , α , (1) which is subject to the following constraints: ∑ i = 1 m r i , j , g = 1 , ∀ g , j , (2) 0 ≤ r i , j , g − ≤ r i , j , g ≤ r i , j , g + < 1 , ∀ g , i , j , (3) | α | 1 = ∑ d = 1 n | α d | ≤ C , (4) where r i , j , g − , r i , j , g + are user-defined intervals for mislabeling risks (based, e.g., on some prior knowledge) and C is the a priori unknown constant that implicitly confines the number of non-zero components of the parameter vector α. The user-defined choice of matrices r ^+/− is not really arbitrary and should be done based on prior knowledge in such a way that the r-constraints (3) do not lead to an empty set.

It is straightforward to verify that the particular case of problem (1–4) with fixed r being an identity matrix is equivalent to the widely used Lasso (or l1 − ) regularization methods introduced by Tibshirani (1996), Bair and Tibshirani (2004), Friedman et al. (2010), Simon et al. (2011), and Taylor and Tibshirani (2015). In context of these Lasso-regularized methods, decreasing the constant C in (4) one reduces the number of non-zero elements in α, thereby reducing the number of non-zero parameters and avoiding overfitting. This becomes especially important in biomedical applications, where the number n of model parameters is large compared with the size of the available statistics—a typical scenario when the danger of overfitting becomes imminent.

As proven in Section 2 of the SI, given the imperfectly labeled datasets X , Y obs , the solution of this optimization problem (1–4) results in the parameter vector α* for any fixed combination of r and C. This solution will be optimal in the case of the log-likelihood, that is, choosing this particular α* will result in a maximal probability for observing the given data X , Y obs , confined to a particular choice of r and C. Selection of the optimal parametric model class ϕ, as well as selection of parameters r and C, can be approached with the standard model selection procedures of machine learning, for example, by means of the cross-validation, with the help of the information criteria or through selective inference (Burnham and Anderson, 2002; Zhang et al., 2010; Taylor and Tibshirani, 2015) (see, e.g., Figure 1A). Uncertainty of the obtained mislabeling risks r* and model parameters α* can be obtained using the common non-parametric bootstrap sampling procedure (Efron and Tibshirani, 1993) (see, e.g., Figures 1C,E, 2A).

If m and N _g are not too large, one can use the following numerical scheme: 1) first, the rectangular domain spanning a set of admissible values for mislabeling matrix elements in (3) and admissible values of C in (4) is sampled (e.g., by means of a uniform equidistant grid), and 2) for every particular grid point (r _s , C _s ), one deploys some standard gradient-based optimization method to solve (1–4). In the second step (ii), one can use an interior-point method or the sequential quadratic programming (Nocedal and Wright, 2006)), performing constrained concave optimization of (2) subject to a constraint (4) only, with fixed values of r _s and C _s . For example, when m = 2 and N _g = 1 (the case emerging in examples 1 and 2), there will be only two independent parameters in r. Together with the scalar dimension for the regularization constant C, this will result in a 3D grid (r _s , C _s ).

With respect to the model function ϕ, optimization of the concave problem (1–4) would only require the evaluation and communication of the function values and the gradients of ϕ with respect to α. This means that a solution of the overall problem (1–4) can be easily integrated into the common software packages for labeled data analysis. Moreover, solutions of problem (1–4) for different particular choices of (r _s , C _s ) can be found completely independent of each other, herewith allowing for a highly scalable (“embarrassingly parallel”) implementation on high-performance computing facilities.

Although the raw datasets used in many types of biomedical data analysis are very large, making sense of them requires statistical methods able to handle data problems where the dimensionality n of their feature spaces is typically orders of magnitude larger than the number T of individuals in the groups. According to central limit theorems, the uncertainty of typical parameter estimation procedures will reduce as T ^1/2 with the growing statistics size T only if the individual data instances in the statistics can be assumed to be independent. The more that dependence between the data instances, the slower will be the rate of this uncertainty reduction with T. In this respect, genomic SNP data pose special problems as they contain a huge fraction of dependencies between SNP pairs that are in a linkage equilibrium with each other. Another source of bias is introduced through the impact of latent/unobserved factors (e.g., in the form of stratification effects and latent variables) and various forms of mislabeling. Mislabeling can bias the results obtained by standard relation measures that dwell on an exact labeling assumption (e.g., it can bias the results of the t-test, chi-square test, Fisher’s exact test, odds ratio, and many other methods). Common methods to overcome this problem, including outlier detection methods and HMMs, either rely on the presence of some subsets that are exactly labeled or introduce additional variables that estimate the probability of every particular data point being an outlier. This significantly increases the dimension of the parameter space and the risk of overfitting. These problems remain an issue even when the numbers of individuals sequenced by platforms such as 23andme approach millions. All these issues (n > > T, violation of the independence assumption, latent impacts, and mislabeling) lead to uncertainty in estimating parameters and make biomedical GWAS applications very challenging. These issues also limit the applicability of advanced Big Data tools such as artificial neuronal networks to problems of this type. A promising direction toward solving these problems can be found in methods based on Lasso regularization ideas first introduced by R. Tibshirani and coworkers (Tibshirani, 1996; Bair and Tibshirani, 2004; Friedman et al., 2010; Simon et al., 2011; Taylor and Tibshirani, 2015), through a robust and computationally efficient shrinkage of the feature space and zeroing out of less relevant feature components. As demonstrated before, introducing a measurement of the probability of group-specific mislabeling and deploying Bayesian tools permit a natural extension of these ideas to situations in which none of the data groups is presumed perfectly labeled. This is accomplished without introducing many new parameters that have to be estimated. For example, in the case of a one-data group (N _g = 1) with two data labels (m = 2, e.g., “Wellderly” and “non-Wellderly” labels in example 2), only two additional mislabeling parameters need to be estimated. The open-source MATLAB implementation we provide here permits implementation of the algorithm in a strongly scalable way. A full analysis of the data in example 2 takes 24 days on a single-core PC, 2 days on a PC workstation with 12 cores, and only 5 h on a small-scale computer cluster with hundred nodes. Measuring the performance of such methods has a general problem due to the presence of latent impacts and mislabeling, which can bias either standard measures such as AUC, accuracy scores, and Fisher’s exact test or and linkage disequilibrium measures. In future studies, a better understanding of ever-growing sets of biomedical data requires the further development of robust and computationally scalable relation measures that can explicitly infer and take into account eventual latent effects and mislabeling.