Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework

Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework View ORCID Profile Bryant Han , Qingling Duan , Ting Hu doi: https://doi.org/10.1101/2025.07.18.664723 Bryant Han 1 Department of Biomedical and Molecular Sciences, Queen’s University , 18 Stuart St, K7L 3N6, ON, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Bryant Han For correspondence: 19bh19{at}queensu.ca Qingling Duan 1 Department of Biomedical and Molecular Sciences, Queen’s University , 18 Stuart St, K7L 3N6, ON, Canada 2 School of Computing, Queen’s University , 25 Union St, K7L 2N8, ON, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ting Hu 2 School of Computing, Queen’s University , 25 Union St, K7L 2N8, ON, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Preview PDF Abstract Motivation Machine learning models in biomedicine have become increasingly complex, often functioning as black boxes. However, understanding contributors to disease and making actionable health interventions requires interpretable models. Common explainable AI methods like SHAP focus on feature importance but fall short in explaining why features contribute in certain patterns or what interventions to take. Counterfactual explanations address this by proposing “what if” scenarios but current tools focus on individual predictions and fail to generalize complex trends. Results We propose the framework Counterfactual Optimization for Actionable interpretabiLity in AI (COALA). COALA interprets models by identifying optimal counterfactuals across user-defined mutable feature subsets and constraining remaining features to reveal how constraint features determine what interventions are optimal. By analyzing counterfactual profiles of features rather than individual features, COALA reveals holistic patterns. Using synthetic and real datasets, COALA reveals simple and complex model trends and provides more intuitive, multi-feature interventions than SHAP. Availability and Implementation Code for COALA implementation, synthetic data, models trained on synthetic data, and code to replicate results and figures are available at https://github.com/brt-solo/COALA . Introduction Many diseases and phenotypes of interest in biomedicine are complex, where multiple genetic factors, environmental factors, or both, impact the outcome ( Virolainen et al., 2023 ). The effects of these factors on complex diseases are not only additive, meaning their effects are independent of each other, but can interact with each other, where the effect one factor has on disease may depend on the presence of a different factor. It is also possible for interactions to occur between more than just two variables and non-linearly ( Motsinger-Reif et al., 2024 ). The biological system of diseases has generally been studied by breaking the system into smaller components and analyzing the parts one at a time, such as by performing multiple association tests between each variable of interest and the outcome ( Motsinger-Reif et al., 2024 ; Virolainen et al., 2023 ). However, knowing that biological systems can be very complex, particularly in understanding diseases and phenotypes, association tests that analyze one or two variables are not sufficient alone to fully understand these diseases. It is necessary to be able to take many variables into account simultaneously, analyzing the system, which can be done through machine learning (ML) ( Xu and Jackson, 2019 ). Explainable artificial intelligence ML in recent years has become increasingly utilized in biomedical research to analyze these complex diseases, while also incorporating high-dimensional datasets ( Xu and Jackson, 2019 ). Although ML models have the potential to capture the complexity of biological patterns, ML models become more uninterpretable with complexity ( Tanyel et al., 2023 ; A. and R., 2023). The interpretation of complex ML models remains a key challenge and barrier in fully embracing artificial intelligence (AI) to make interpretations and new hypotheses. Even if an ML model can more accurately predict a disease outcome, the disease cannot be better understood without first understanding the model. This necessity to understand ML models across many applications has given rise to the field of explainable AI (XAI) (A. and R., 2023). The need for XAI in biomedicine is especially apparent, as healthcare applications account for the largest proportion of XAI-related publications (A. and R., 2023). XAI is a rapidly developing field with diverse approaches that usually fall into two general categories: model-focused and output-focused interpretability (A. and R., 2023). Model- focused approaches aim to interpret the mechanisms of a complex ML model, such as in neural networks ( Bereska and Gavves, 2024 ). This is analogous to understanding the structures in a brain and the role each lobe or neuron has. Output-focused approaches, on the other hand, focus on explaining the predictions by an ML model (A. and R., 2023). Output-focused approaches can be categorized into intrinsically interpretable models and post-hoc explanations (A. and R., 2023). Intrinsically interpretable models are self-explanatory and have a simple structure, such that no additional analysis is required to fully understand it. However, such simpler models are often not sufficient for accurately predicting complex outcomes. Hence, many post-hoc explanation tools have been developed, which are also widely used for ML models in biomedicine (A. and R., 2023; Chamola et al., 2023 ). The dominant post-hoc methods in biomedical literature are attribution based, such as Locally Interpretable Model- Agnostic Explainer (LIME) and Shapley Additive exPlanations (SHAP), which assign scores to each feature’s contribution to an individual prediction ( Ribeiro et al., 2016 ; Lundberg and Lee, 2017 ). The scores from feature attribution can then be aggregated to produce a feature importance score for each feature ( Molnar, 2025 ). While feature attribution can explain what contributes to predictions, it does not necessarily explain why that feature contributes to such degree for that prediction, especially if there are interactions between variables and the relationships between variables or between variables and outcome are non-linear ( Molnar, 2025 ). Feature attribution is suitable for explaining the present, but not for informing on interventions in the future. One post-hoc method that can address the limitations of feature attribution is through counterfactuals ( Byrne, 2019 ). Counterfactuals are hypothetical instances modified from original instances that would lead to a different outcome. A counterfactual statement would follow “If X variable is changed by Y, the outcome changes by Z”. This creates a direct link between a variable and the outcome, clearly showing that X variable causes a change in the predicted outcome. While the inferred causality is about the model, and cannot be used to define real-world causality, ML models may reveal novel patterns that can inform on real patterns ( Baron, 2023 ; Molnar, 2025 ). Counterfactual approaches are particularly promising not just because they reveal causality within the model, but also because they draw directly from human psychology. As counterfactuals are a common component of human reasoning, explanations about an ML model’s predictions can be intuitive to non-AI experts in biomedicine ( Byrne, 2019 ; Nagesh et al., 2023 ). However, current counterfactual approaches are not without their own limitations. First, counterfactual tools have generally been applied to explaining individual instances ( Wielopolski et al., 2024 ). This becomes difficult to rely on if the ML model has a relatively low accuracy, which is often the case in biomedicine ( Papenmeier et al., 2019 ). Even if an ML model is more accurate than human predictions, it cannot be certain a specific counterfactual is trustworthy. Additionally, counterfactuals of specific instances cannot necessarily be extrapolated to all instances, making a counterfactual’s insight very limiting ( Molnar, 2025 ). Thus, to overcome this limitation, counterfactuals should be used not just to explain specific instances, but to explain trends on a population level ( Wielopolski et al., 2024 ; Brughmans et al., 2024 ). Second, variables can be interconnected in how they affect complex diseases, or more accurately, how they affect the prediction ( Virolainen et al., 2023 ). Thus, counterfactuals should be systematically generated by changing multiple features at a time to understand how multiple features cumulatively change the prediction outcome. Third, for a given instance, there are multiple counterfactuals that can result in the same change of outcome. As such, there is often a seemingly limitless number of possible counterfactuals ( Brughmans et al., 2024 ). Thus, there needs to be an objective of optimizing counterfactuals, such as maximizing the predicted outcome. Existing counterfactual approaches, such as DiCE ( Mothilal et al., 2020 ) and MACE ( Yang et al., 2022 ), are designed to generate plausible counterfactuals for individual instances rather than to analyze patterns between populations. As such, they do not support subgroup identification or explain how optimal interventions vary across individuals. In this research, we propose the first counterfactual framework explicitly developed to provide population-level interpretability by identifying how optimal counterfactuals are modulated by constraint features across a dataset. Evolutionary algorithms To optimize for counterfactuals that lead to the best outcomes, we use evolutionary algorithms, which are search algorithms inspired by biological evolution and natural selection ( Yu and Gen, 2010 ). Evolutionary algorithms are suitable for optimizing problems that have similar challenges to nature. For problems with a small search space, a random search or grid search may be sufficient and problems that are differentiable and have a single global optima may be solved with gradient descent ( Yu and Gen, 2010 ). However, in nature, populations can be very diverse, and it is possible for high fitness individuals to be very different from each other yet similar to some low fitness individuals, just like many optimization problems ( Sudholt, 2018 ). Evolutionary algorithms have three main characteristics that operate similarly to natural evolution: population, fitness, and variation ( Yu and Gen, 2010 ; Lutton et al., 2016 ). First, there is a population of individuals. For our problem, each individual, which is a candidate solution, is a generated counterfactual and a population of candidates is generated in the process of finding optimal counterfactuals for a single sample. Fitness, which is the metric for selection, is defined as the value of the ML model-predicted outcome. Variation’s purpose is to generate new candidate solutions that may have a higher fitness, which can be done most commonly through mutation (random permutation of a candidate) and crossover (combining two parent candidates to produce an offspring candidate). A successful approach to real-life challenges is creating diversity. Whether it be a species attempting to survive in a crowded habitat or an exam that tests a student’s problem solving skills, only trying to move towards what seems like the best solution at the moment is not always the best, nor the fastest, as other high fitness solutions may be missed ( Sudholt, 2018 ; Markert et al., 2010 ). The Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) algorithm is an evolutionary algorithm that uses diversity to illuminate the search space ( Mouret and Clune, 2015 ). This is achieved by partitioning the feature space of possible solutions into a grid of cells, for which MAP-Elites identifies the optimal solution for each cell. Counterfactuals in biomedicine Counterfactuals are used constantly in medicine, such as when a physician prescribes a treatment or suggests a lifestyle change in order to improve a health outcome (Béal and Latouche, 2020 ). Physicians may also aim to provide personalized treatments, as not all patients will benefit equally from the same recommendations. The ideal treatment for a patient depends on their other health-related variables ( Virolainen et al., 2023 ; Hamburg and Collins, 2010 ). Regarding our counterfactual approach, changes made for an optimal counterfactual depend on the values of the unchanged variables. Only for simple diseases would all individuals benefit equally from identical treatments. Thus, we shape our counterfactual analysis framework around this aspect of complex diseases. Other XAI tools, such as SHAP or other counterfactual tools, have suggested how a feature of interest affects the outcome, but we focus on understanding the interactions and dependencies between features ( Lundberg and Lee, 2017 ; Mothilal et al., 2020 ). Framework overview Here we introduce Counterfactual Optimization for Actionable interpretabiLity in AI (COALA), a novel framework that, given a population of samples, searches for optimal counterfactuals to provide interpretable explanations of the model and potential actionable interventions ( Figure 1 ). COALA draws from the MAP-Elites algorithm to generate a diverse set of optimal counterfactuals ( Mouret and Clune, 2015 ). COALA partitions the search space into cells that are defined by what features were mutable to generate the counterfactual and finds the optimal counterfactual in each cell. Download figure Open in new tab Fig. 1. Overview of the pipeline for analyzing counterfactuals with COALA. Values of specific feature categories are changed for an input vector, with high fitness candidates placed in its corresponding cell, where the fitness is the ML model-predicted outcome. Optimal counterfactuals are expected to vary if there are interactions between constraint and mutable features, such that identified clusters represent groups of patients and their appropriate intervention. COALA offers several advantages over existing XAI and counterfactual methods. Unlike instance-specific approaches, COALA focuses on population-level trends, revealing how patterns generalize across individuals. By systematically changing groups of features, it supports multi-feature counterfactual analysis, uncovering complex but interpretable interventions. Additionally, by identifying the optimal counterfactuals in each group, COALA facilitates direct comparisons across individuals and highlights how the most effective interventions differ between subpopulations. Methods Counterfactual optimization framework COALA is conceptually simple and intuitive to implement and extend. The user first supplies a pre-trained machine learning model for a prediction task (e.g., disease risk), which is used to evaluate candidate solutions based on their predicted outcomes. Next, the user defines feature categories (lines 2–4 of Algorithm 1 ), which serve to partition the search space into interpretable cells. Each cell then represents a combination of feature pairs that can be changed together to create a counterfactual. These feature categories can be specified based on domain-specific knowledge. For instance, in a model trained on electronic health records (EHR), features might be grouped into categories such as diagnoses, demographics, and medications. In a multi-omics context, each omic layer (e.g., transcriptomics, proteomics, metabolomics) could form a separate category. A candidate is defined as a generated counterfactual, and the fitness of that candidate is defined as the value of the ML model-predicted outcome. An optimal counterfactual is the candidate with the highest fitness. For each input sample x , COALA algorithmically identifies an optimal counterfactual for each cell of the feature map. Using the input sample as a fixed reference point, COALA begins with a random initialization phase (lines 7–10 and 16–20 of Algorithm 1 ), where a diverse population of counterfactuals is generated. In each iteration of this phase, a cell ( 𝒞 i,j ) is selected at random, and a counterfactual is generated by randomly perturbing the features in the corresponding category pair ℱ i,j . In our basic implementation of COALA, each continuous mutable feature x k in the feature set ℱ i,j is randomly sampled from a uniform distribution bounded within three standard deviations of the mean, ensuring that all counterfactuals remain within the observed range: Binary features are randomly assigned values from their respective two-category sets. Repeating this process produces an initial population of candidate solutions across cells. During the optimization phase (lines 7–8 and 12–20 of Algorithm 1 ), COALA continues to iteratively refine counterfactuals. In each iteration, a cell is selected at random, and two counterfactuals from that cell—one elite solution ( 𝒳 i,j ) and one sampled at random—are selected as parents. These are recombined via crossover to produce a new candidate, which is assigned back to the same cell, regardless of whether it outperforms its parents. This promotes local exploration within the cell’s defined feature subspace. For our analyses, we used uniform crossover, but COALA also supports single point crossover, simulated binary crossover, and random mutation. The optimization proceeds until a termination criterion is met, such as a maximum number of generations. This process is repeated independently for each input sample. Countless combinations of feature groups can be used to define multiple cells, but for simplicity, we first used COALA on the models and datasets by limiting the counterfactual search space to one cell, where only real-life fixed variables (e.g., genetic predisposition) are constrained and real-life actionable features (e.g., nutrition) are mutable. After using COALA by only changing actionable features, we then used domain knowledge to define multiple feature categories and cells for the synthetic data and real dataset for further COALA analysis. Actionable interpretations For each sample, COALA will generate one counterfactual for each of I iterations. Each sample will have one optimal counterfactual for each cell, so if the search space was split into ten cells, there would be ten optimal counterfactuals. Each cell’s counterfactuals can then be analyzed with common statistical methods to identify trends in counterfactuals across the sample population. This involves clustering of the counterfactuals to identify different risk groups, then statistical tests such as ANOVA to identify risk factors that drive a sample’s risk group membership. Because COALA searches for optimal counterfactuals within specific cells, where all features in ℱ i,j are mutable and all other features are held constant, differences in the resulting optimal counterfactuals across samples can be attributed to differences in the constraint features. In other words, COALA identifies the optimal values of mutable features conditional on the fixed values of the constraint features. If the relationship between features and the outcome is simple and additive, one would expect similar mutable feature changes to be optimal across all samples. Conversely, if counterfactuals differ across samples, this suggests that the model has captured interactions between mutable and constraint features. To identify risk groups, samples are clustered based on the mutable features of their optimal counterfactuals within a given cell. This reveals distinct clusters of counterfactuals, which may be candidate interventions, and the clustering process is done independently once for each cell. Although various clustering methods may be appropriate, we performed K-means clustering and then principal component analysis (PCA) for visualization. The optimal number of clusters was selected using the silhouette score, which quantifies the cohesion and separation of clusters. Algorithm 1 COALA: Constraint-Aware Optimization for Local Adjustments Download figure Open in new tab We then identified which features characterize the different treatment clusters by performing ANOVA tests on the mutable features and computing the proportion of each feature’s variance explained by cluster membership ( η 2 ): We also performed ANOVA tests on the constraint features of each cell and compute η 2 . All ANOVA tests were conducted on the original (unscaled) feature values to preserve interpretability. Constraint features with larger relative η 2 values can be interpreted as interacting with mutable features and modulating their influence on the predicted outcome. In other words, these constraint features act as drivers in determining optimal treatments. SHAP analysis We compared COALA with a SHAP-based clustering framework commonly used to identify risk subgroups in machine learning models. SHAP values were computed for each sample, yielding a vector of feature attributions. Samples were then clustered with hierarchal clustering using Ward’s linkage method and Euclidean distance based on their SHAP vectors to identify subgroups with similar patterns of model explanation. The number of clusters was specified based on the expected number of clusters from COALA analysis. We took our SHAP clustering approach from biomedical research that use ML models to identify actionable interventions for complex diseases( Chen et al., 2025 ). To interpret the resulting clusters, ANOVA tests were performed on both the SHAP values and original feature values to identify which features differentiated the identified subgroups the most. Datasets and models We compared COALA and SHAP-based clustering on a real biomedical dataset and a synthetic dataset, both of which are tabular datasets with a continuous outcome. Due to restrictions and privacy concerns of the biomedical dataset, exact feature names, outcome, and details about the cohort cannot be disclosed. The features were categorized prior to using COALA based on domain knowledge of biological systems. To benchmark the COALA framework in a controlled setting, we constructed a synthetic dataset with n = 1000 samples, 9 input features divided into four conceptual views: genetic risk, environmental exposures, nutritional, and metabolic, and a continuous outcome y that represents a hypothetical health measurement (e.g., cholesterol levels) ( Table 1 ). View this table: View inline View popup Download powerpoint Table 1. Structure of the synthetic dataset by feature view. The outcome variable y was generated using a weighted combination of additive effects and interaction terms: A linear regression model and a multilayer perceptron (MLP) were each trained on the synthetic dataset and an MLP was trained on the real dataset, using an 80/20 train–test split. Linear regression models are limited to capturing additive relationships, while MLPs are capable of modeling nonlinearities and feature interactions. Each trained model–dataset pair was then analyzed using both the COALA and SHAP frameworks. We also analyzed a ground truth model created to perfectly match the synthetic data to assess whether COALA can identify true patterns. Results Evaluation on synthetic datasets We first evaluated COALA on a linear regression model and an MLP, both trained on the same synthetic dataset, to investigate if COALA can identify the same optimal counterfactual for each sample using different models. As expected, COALA was able to identify identical optimal changes for all samples when applied to the linear model. When applied to the ground truth model, there were 8 distinct optimal counterfactuals ( Figure 2A ). ANOVA test of the counterfactuals’ mutable feature values revealed that the variance in mutable features E1, M2, and N2 can be explained by cluster membership ( Figure 2B ), and the variance in mutable features can be explained by the variance in the constraint feature values G1 and G2 ( Figure 2C ). In other words, COALA analysis of the ground truth model showed that G1 and G2 dictate what the optimal values of E1, M2, and N2 should be. Additionally, as expected E2, N1, M1, and G3 have no variation across samples as they were not specified to have any interactions during data generation 3. Download figure Open in new tab Fig. 2. COALA analysis of the ground truth model for the synthetic dataset. (A) PCA of mutable features across counterfactuals colored by cluster. (B) Proportion of variance in mutable features explained by cluster membership ( η 2 ). (C) Proportion of variance in constraint features explained by cluster membership ( η 2 ). (D) Mean values of mutable features by cluster. (E) Mean values of constraint features by cluster. By examining the mean values of features in each cluster, we observed that the counterfactual for E1 is positive when the constraint G1 is positive, and negative when G1 is negative ( Figure 2D, E ). This is expected from the specified G1 · E1 interaction during data generation ( Equation 3 ). We also observed that when G2 was positive, N2 and M2 were either both positive, or both negative. When G2 was negative, one of N2 and M2 was positive, and the other negative. This pattern matches the expected results from the G2 · N2 · M2 interaction. The synthetic dataset was split into an 80/20 split for training and validation of the MLP model, and standardized prior to training. The MLP model had an MSE of 0.4289 and R 2 of 0.9609. There were seven distinct clusters identified ( Figure 3A ), with N2, M2, and E1 and G2 and G1 having the largest effect sizes among mutables and constraints, respectively ( Figure 3B, C ). This matches the COALA analyses of the ground truth model. However, E2 is also shown to have a large effect size, which is not an expected trend ( Figure 3B ). Although E2 has no interactions with other variables, the MLP model overfits and identifies Cluster 4 to benefit from a lack of E2 ( Figure 3D ). Thus, COALA can also be used to confirm if a model is correctly identifying trends. Download figure Open in new tab Fig. 3. COALA analysis of the MLP model trained on the synthetic dataset. (A) PCA of mutable features across counterfactuals colored by cluster. (B) Proportion of variance in mutable features explained by cluster membership ( η 2 ). (C) Proportion of variance in constraint features explained by cluster membership ( η 2 ). (D) Mean values of mutable features by cluster. (E) Mean values of constraint features by cluster. Application to a real dataset The MLP trained on the real biomedical dataset had an R 2 of 0.0860 when applied to the validation set. Although it a has lower accuracy than clinically used ML models, COALA still reveals consistent learned trends even if the final prediction is not the most precise. The silhouette score identified three clusters as the optimal clustering of counterfactuals for the cell where only real-life fixed features are constrained, although visually PCA does no reveal distinctly compact grouping ( Figure 4A ). However, ANOVA test identified features Actionable 1, 3, 9, 2, and 8 to vary across clusters, with Constraint 1 being the main driver of these differences ( Figure 4B, C ). Cluster 1 counterfactuals notably have a mean positive Actionable 1 and 3, whereas only the mean of Actionable 3 is positive for Cluster 2 and Actionable 1 is positive for Cluster 3 ( Figure 4D ). Results suggest that the different optimal counterfactuals are driven by Fixed 1, which is low in Cluster 2 and high in Cluster 3 ( Figure 4C, E ). These results are supported by the previously established trend that Fixed 1 forms interactions with the actionable features and is a key risk factor. ML models used to predict complex diseases can suffer from low accuracy, but this also shows that COALA can be used to confirm if a model recognizes expected trends. Download figure Open in new tab Fig. 4. COALA analysis of the MLP model trained on a real biomedical dataset. (A) PCA of mutable features across counterfactuals colored by cluster. (B) Proportion of variance in mutable features explained by cluster membership ( η 2 ) The red dotted line indicates the η 2 threshold for a statistically significant ANOVA test (p¡0.05). (C) Proportion of variance in constraint features explained by cluster membership ( η 2 ). (D) Mean values of mutable features by cluster. (E) Mean values of constraint features by cluster. Because clustering of the counterfactuals was not distinct, it suggests that there are not discrete optimal counterfactuals, or treatments, but rather in this case fall on a continuous spectrum. This is supported by the variance in counterfactuals following the gradient of Constraint 1 values ( Figure 5 ). It can be seen that Fixed 1, though important, is not the sole driver of optimal mutable features, as some samples with relatively higher Fixed 1 can be seen to have more similar counterfactuals with low Fixed 1 samples and vice versa. This leaves much to explore with how COALA can be used to extract more detailed models of how features interact. Download figure Open in new tab Fig. 5. COALA analysis of the MLP model trained on a real biomedical dataset visualized with PCA of mutable features across counterfactuals, colored by the value of Constraint 1. Multiple cells We applied COALA to the ground truth model of the synthetic dataset and MLP model trained on a real dataset using multiple cells ( Figure 6 , 7 ). This allowed for the identification of specific feature pairs that interact with each other in the synthetic dataset, such as identify E1 as a mutable that is driven by G1 when only Environment features are mutable ( Figure 6B, C ). There were cells where no constraint was identified to reveal the G2 · N2 · M2 interaction ( Figure 6C ), because although there are two clusters of optimal counterfactuals for M2 ( Figure 6B ), different groups can share the same counterfactual (e.g., negative G2 and N2 benefit from the same change as positive G2 and N2). However, a simple decision tree that predicts cluster label from constraint values can untangle the different constraint profiles present in each cluster ( Figure 8A ). Download figure Open in new tab Fig. 6. COALA analysis of ground truth model for the synthetic dataset, with multiple cells specified. (A) PCA of mutable features across counterfactuals colored by cluster. (B) Proportion of variance in mutable features explained by cluster membership ( η 2 ). (C) Proportion of variance in constraint features explained by cluster membership ( η 2 ). Download figure Open in new tab Fig. 7. COALA analysis of the MLP model trained on the biomedical dataset, with multiple cells specified. (A) PCA of mutable features across counterfactuals colored by cluster. (B) Proportion of variance in mutable features explained by cluster membership ( η 2 ). The red dotted line indicates the η 2 threshold for a statistically significant ANOVA test (p¡0.05). (C) Proportion of variance in constraint features explained by cluster membership ( η 2 ). Download figure Open in new tab Fig. 8. Decision trees predicting cluster membership from constraint features for one cell of two COALA experiments. (A) Ground truth model for the synthetic dataset, where only metabolic features are mutable. (B) MLP model trained on a real dataset, where only Actionable Category 2 features are mutable. When only features of Actionable Category 2 were mutable in the analysis of the MLP trained on a real dataset, features Fixed 1, Actionable 5, and Actionable 6, were identified as leading drivers for the optimal Mutable 9 and 10 ( Figure 7B, C ). When all Actionable features are mutable, Fixed 1 was identified as the clear main driver, but Actionable features may be just as if not more important than Fixed 1 if additional features are constrained. It can be seen across cells in Figure 7C how the Constraint 1 varies in its relative effect size compared to other constraints across cells, emphasizing the importance of optimizing for multiple cells. How the constraints collectively affect what the optimal counterfactual, however, is not as simple as the synthetic dataset ground truth, as a decision tree is too simple to consistently separate clusters based on the constraints ( Figure 8B ). The cell for when only Actionable Category 2 is mutable reveals only 4 distinct solutions, so it can be expected that samples can be relatively easily separated into clusters based on the constraints, such as with a decision tree. However, because this is not the case, this suggests that even if two samples have similar or identical optimal counterfactuals, they do not necessarily have similar constraints. The opposite may also be true, where two samples have contrasting counterfactuals, yet relatively similar constraints. Comparison with SHAP Using SHAP, Fixed 1 was identified as having the greatest impact on the predicted outcome. Although it may be hypothesized that Fixed 1 forms interactions with other features, what those features may be cannot be inferred from SHAP analysis alone. Actionable 8 and 4 were also identified to have a high feature importance. We applied SHAP clustering to the MLP model trained on a real biomedical dataset. While neither COALA nor SHAP clustering reveal distinct risk groups ( Figure 4A , 9A ), SHAP clustering fails to identify Fixed 1 as an important feature when determining an optimal treatment ( Figure 9B, C ). Although Fixed 1 has a high feature importance, the impact it has on the predicted outcome does not vary enough across samples for it to be important in clustering. Using SHAP clustering identified Actionable 10, 8, and 9 as features that define risk groups ( Figure 9B, D ). The variance in SHAP values among features generally mirrors the variance in their original input values ( Figure 9C, E ). This suggests that the differences in the features’ SHAP values were primarily due to the difference in each individual feature’s original values, and not by setting a context that determines the effects from other features. SHAP is still informative, showing which features have impacted different predictions, but still leaves questions of why and how( Figure 9D, E ). SHAP alone was not able to show if a feature’s effects when changed would vary across samples, whereas COALA showed how the constraints set the context of how the mutable features’ impact the outcome. Download figure Open in new tab Fig. 9. SHAP clustering analysis of the MLP model trained on a real biomedical dataset. (A) PCA for visualization of samples’ SHAP values. (B) Proportion of variance in SHAP values explained by cluster membership ( η 2 . (C) Proportion of variance in original feature values explained by cluster membership ( η 2 ). (D) Mean SHAP values of feature by cluster. (E) Mean original values of feature by cluster. Discussion The main advantage of COALA is how it frames features into constraints that act as drivers on the mutable features. How COALA identifies features of interest is not only intuitive, but how it identifies risk subgroups is more closely aligned to clinical interests and personalized intervention strategies. On the other hand, SHAP-based clustering identifies groups of individuals who share similar reasons for their prediction, grouping together samples with similar feature attributions. For instance, individuals in one cluster may be strongly influenced by a specific feature (e.g., Actionable 1 and Fixed 1, Figure 4B, C ). However, SHAP could not reveal if individuals outside that cluster would respond similarly to that same feature. SHAP identified risk groups based on features that directly contribute to the predicted outcome, clustering individuals who share similar influential features ( Figure 9B, C ). In contrast, COALA identified risk groups by examining how certain features modulate the influence of other features on the outcome, revealing interaction effects ( Figure 4B, C ). In complex systems, outcomes are often shaped by interactions between features, where the impact of one variable depends on the context set by others. SHAP assesses an instance retrospectively, and while it may give insights to what changes could be made, how beneficial it would be to improving the predicted outcome is not always clear. COALA is designed to assess instances prospectively, and provides more certainty to how changes would impact the predicted outcome. There are pre-existing interaction models that are commonly used, but they assess pairwise effects independently, limiting the potential complexity that can be captured ( Motsinger-Reif et al., 2024 ). ML models have the potential to capture such complexity, but interpretations will similarly remain limited if XAI methods analyze one or two features at a time. Thus, COALA compliments the strength of complex ML models, that is to use high-dimensional inputs, by accounting for multiple interacting features, enabling it to capture higher- order interactions. Although COALA was designed and applied to ML models trained on tabular data, COALA’s approach of using counterfactuals is not limited to any specific type of dataset, such as EHR data and multiomic data. COALA is also model agnostic, as it treats all models as a black box. It only needs the prediction function to optimize for counterfactuals. For interpretable models, such as linear regression and K-nearest neighbour (KNN) models, that already provide transparent explanations for users to generate counterfactuals on their own, COALA is not necessary. However, application to black box models, such as MLPs in this paper, is when COALA can thrive. COALA is able to map and differentiate samples based on what the optimal counterfactual is, but as expected, how cluster membership can be determined with the constraints is not necessarily simple ( Figure 8 ). Future applications can explore nuanced approaches to understand how constraints collectively determine what the optimal mutables are. Counterfactuals by nature imply causality, but it is important to note that in XAI the causality directly applies only to the model, and not to the real-world relationships, as such causal inferences require causal models ( Molnar, 2025 ). However, causal models often struggle with more complex patterns, such as nonlinear interactions and context-specific effects (Runge et al.). In general, causal modeling requires strong domain knowledge of the features and outcome to make assumptions for causal modeling, which can be quite challenging with a high dimensional dataset for a complex outcome ( Berkessa et al., 2025 ). Even though COALA can not fully validate if trends identified by an ML model are truly causal, it can help make new hypotheses and inform users of ways to build upon a causal model. COALA’s objective is not to replace established interaction analysis methods or causal models, but to address the necessity of methods that compliment black box models. Author contributions statement B.H., Q.D., and T.H. conceived the experiments, B.H. conducted the experiments, analyzed the results, and wrote the manuscript, B.H., Q.D., and T.H. reviewed the manuscript. Funder Information Declared Canadian Institutes of Health Research, https://ror.org/01gavpb45 References S. A. and S. R. A systematic review of explainable artificial intelligence models and applications: Recent developments and future trends . Decision Analytics Journal , 7 : 100230 , 2023 . doi: 10.1016/j.dajour.2023.100230 . OpenUrl CrossRef PubMed ↵ S. Baron . Explainable AI and causal understanding: Counterfactual approaches considered . Minds and Machines , 33 ( 2 ): 347 – 377 , 2023 . doi: 10.1007/s11023-023-09637-x . OpenUrl CrossRef ↵ L. Bereska and E. Gavves . Mechanistic interpretability for AI safety – a review , 2024 . ↵ Z. A. Berkessa , E. Läärä , and P. Waldmann . A review of causal methods for high-dimensional data . IEEE Access , 13 : 11892 – 11917 , 2025 . doi: 10.1109/ACCESS.2024.3524261 . OpenUrl CrossRef ↵ D. Brughmans , L. Melis , and D. Martens . Disagreement amongst counterfactual explanations: how transparency can be misleading . TOP , 32 ( 3 ): 429 – 462 , 2024 . doi: 10.1007/s11750-024-00670-2 . OpenUrl CrossRef ↵ R. M. J. Byrne . Counterfactuals in explainable artificial intelligence (XAI): Evidence from human reasoning . In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence , pages 6276 – 6282 , Macao, China , 2019 . International Joint Conferences on Artificial Intelligence Organization . ISBN 978-0-9992411-4-1 . doi: 10.24963/ijcai.2019/876 . OpenUrl CrossRef ↵ J. Béal and A. Latouche . Causal inference with multiple versions of treatment and application to personalized medicine , 2020 . arxiv: 2005.12427 [stat]. ↵ V. Chamola , V. Hassija , A. R. Sulthana , D. Ghosh , D. Dhingra , and B. Sikdar . A review of trustworthy and explainable artificial intelligence (XAI) . IEEE Access , 11 : 78994 – 79015 , 2023 . doi: 10.1109/ACCESS.2023.3294569 . OpenUrl CrossRef ↵ S.-F. Chen , S. E. Lee , H. J. Sadaei , J.-B. Park , A. Khattab , J.-F. Chen , C. Henegar , N. E. Wineinger , E. D. Muse , and A. Torkamani . Meta-prediction of coronary artery disease risk . Nature Medicine , pages 1 – 12 , 2025 . doi: 10.1038/s41591-025-03648-0 . Publisher: Nature Publishing Group . OpenUrl CrossRef ↵ M. A. Hamburg and F. S. Collins . The path to personalized medicine . New England Journal of Medicine , 363 ( 4 ): 301 – 304 , 2010 . doi: 10.1056/NEJMp1006304 . OpenUrl CrossRef PubMed Web of Science ↵ S. M. Lundberg and S.-I. Lee . A unified approach to interpreting model predictions . In Advances in Neural Information Processing Systems , volume 30 , 2017 . ↵ E. Lutton , N. Perrot , and A. Tonda . A Brief Introduction to Evolutionary Algorithms, chapter 2 , pages 7 – 32 . John Wiley Sons, Ltd , 2016 . ISBN 9781119136828 . doi: 10.1002/9781119136828.ch2 . OpenUrl CrossRef ↵ J. A. Markert , D. M. Champlin , R. Gutjahr-Gobell , J. S. Grear , A. Kuhn , T. J. McGreevy , A. Roth , M. J. Bagley , and D. E. Nacci . Population genetic diversity and fitness in multiple environments . 10 : 205 , 2010 . ISSN 1471-2148 . doi: 10.1186/1471-2148-10-205 . OpenUrl CrossRef PubMed ↵ C. Molnar . Interpretable Machine Learning . 3 edition, 2025 . ↵ K. Mothilal , A. Sharma , and C. Tan . Explaining machine learning classifiers through diverse counterfactual explanations . In Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency , pages 607 – 617 , Barcelona Spain , 2020 . ACM . ISBN 978-1-4503-6936-7 . doi: 10.1145/3351095.3372850 . OpenUrl CrossRef ↵ A. A. Motsinger-Reif , D. M. Reif , F. S. Akhtari , J. S. House , C. R. Campbell , K. P. Messier , D. C. Fargo , T. A. Bowen , S. S. Nadadur , C. P. Schmitt , K. G. Pettibone , D. M. Balshaw , C. P. Lawler , S. A. Newton , G. W. Collman , A. K. Miller , B. A. Merrick , Y. Cui , B. Anchang , Q. E. Harmon , K. A. McAllister , and R. Woychik . Gene-environment interactions within a precision environmental health framework . Cell Genomics , 4 ( 7 ): 100591 , 2024 . doi: 10.1016/j.xgen.2024.100591 . OpenUrl CrossRef PubMed ↵ J.-B. Mouret and J. Clune . Illuminating search spaces by mapping elites , 2015 . ↵ S. Nagesh , N. Mishra , Y. Naamad , J. M. Rehg , M. A. Shah , and A. Wagner . Explaining a machine learning decision to physicians via counterfactuals , 2023 . ↵ A. Papenmeier , G. Englebienne , and C. Seifert . How model accuracy and explanation fidelity influence user trust in AI . 2019 . ↵ M. T. Ribeiro , S. Singh , and C. Guestrin . “Why should I trust you?” Explaining the predictions of any classifier . In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining , pages 1135 – 1144 . ACM , 2016 . J. Runge , P. Nowack , M. Kretschmer , S. Flaxman , and D. Sejdinovic . Detecting and quantifying causal associations in large nonlinear time series datasets . 5 ( 11 ): eaau4996 . doi: 10.1126/sciadv.aau4996 . Publisher: American Association for the Advancement of Science . OpenUrl FREE Full Text ↵ D. Sudholt . The benefits of population diversity in evolutionary algorithms: A survey of rigorous runtime analyses , 2018 . ↵ T. Tanyel , S. Ayvaz , and B. Keserci . Beyond known reality: Exploiting counterfactual explanations for medical research , 2023 . ↵ S. J. Virolainen , A. VonHandorf , K. C. M. F. Viel , M. T. Weirauch , and L. C. Kottyan . Gene–environment interactions and their impact on human health . Genes & Immunity , 24 ( 1 ): 1 – 11 , 2023 . doi: 10.1038/s41435-022-00192-6 . Publisher: Nature Publishing Group . OpenUrl CrossRef PubMed ↵ P. Wielopolski , O. Furman , J. Stefanowski , and M. Zieba . Unifying perspectives: Plausible counterfactual explanations on global, group-wise, and local levels , 2024 . arxiv: 2405.17642 [cs]. ↵ C. Xu and S. A. Jackson . Machine learning and complex biological data . Genome Biology , 20 : 76 , 2019 . doi: 10.1186/s13059-019-1689-0 . OpenUrl CrossRef PubMed ↵ W. Yang , J. Li , C. Xiong , and S. C. H. Hoi . MACE: An efficient model-agnostic framework for counterfactual explanation , 2022 . ↵ X. Yu and M. Gen . Introduction to Evolutionary Algorithms . Springer , 2010 . View the discussion thread. Back to top Previous Next Posted July 18, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework Bryant Han , Qingling Duan , Ting Hu bioRxiv 2025.07.18.664723; doi: https://doi.org/10.1101/2025.07.18.664723 Share This Article: Copy Citation Tools Identifying intervention strategies from machine learning models with COALA: a counterfactual optimization framework Bryant Han , Qingling Duan , Ting Hu bioRxiv 2025.07.18.664723; doi: https://doi.org/10.1101/2025.07.18.664723 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41911) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13371) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40364) Molecular Biology (17163) Neuroscience (88537) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15129) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)