Modern machine learning methods allow one to leverage complex relationships present in data to generate accurate predictions but do not reveal them to the investigator. We propose an automated conjecturing framework for discovering nonlinear and boolean relationships among the features in a given dataset. In many situations such as medical and financial decision making, knowing the basis for predictions is important for reasons including understanding causal relationships, establishing trust, and assigning liability (Stoyanovich et al, 2020) . The relationships that are discovered with our framework can be used as features in learning algorithms to boost performance while maintaining model transparency.
In situations where black box models provide accurate predictions, there could be easily-understandable feature relationships that guide system behavior. Our goal is to develop a framework for revealing the least complex nonlinear relationships. For example, consider analyzing measurements on the gravitational force between several pairs of masses at various distances apart. Below we simulate such a system and demonstrate how one can discover that force is directly proportional to the product of the masses and inversely proportional to the square of the distance.
Our framework is based on Fajtlowicz's Dalmation Heuristic (Fajtlowicz, 1995) . The heuristic was originally implemented in Graffiti (Fajtlowicz, 1995; Larson and Van Cleemput, 2017) which was the first program to produce research conjectures that led to new mathematical theory. The program produces statements that are relations between mathematical invariants which are numerical attributes of examples. Recent implementations of the Dalmation Heuristic have been applied to the discovery of relationships for graphs (Larson and Van Cleemput, 2016) and game strategies (Bradford et al, to appear) . The heuristic was adapted to work with properties which are boolean attributes of examples (Larson and Van Cleemput, 2017) .
The proposed framework builds on and fills gaps in computational scientific discovery, interpretable artificial intelligence (AI), automated feature engineering, and statistical methodology.
Computational scientific discovery. Schmidt and Lipson (2009) extend work in symbolic regression to develop a system for discovering laws for dynamical systems by considering relationships among derivatives. Their work led to the development of a software Eureqa. Symbolic regression is the use of genetic programming to approximate a target function on training data and generalize to produce predictions on new data (Nicolau and Agapitos, 2020) and until the work of Schmidt and Lipson, the focus was on prediction by approximating an underlying function rather than characterizing true functional relationships. Recently, Udrescu and Tegmark (2020) combined a variety of strategies including dimensional analysis, symmetry identification, neural network training, and brute-force enumeration into a framework called AI Feynman to recover true physical functional forms from data.
Other frameworks have been proposed to computationally generate conjectures from data and discover scientific laws. Data smashing is introduced by Chattopadhyay and Lipson (2014) as a method for computing dissimilarities from streams of data (e.g., electroencephalogram data) to aid in revealing relationships among observations. Jantzen (2020) proposes an algorithm with the similar purpose of detecting types of dynamical systems called dynamical kinds. Subsequently, these kinds "are then targets for law-like generalization" (Jantzen, 2020). While Jantzen's work provides a method for discovering the kinds, it does not suggest how to recover the "laws". It is these relationships that we aim to discover with the proposed conjecturing framework.
Our work is distinguished from these previous works in that 1) we focus on generating bounds for invariants that serve as hypotheses for the investigator rather than recover true functional forms or generate accurate predictions, 2) our invariant conjecturing algorithm is paired with a property conjecturing algorithm for discovering both nonlinear and boolean relationships 3) our framework is designed for a given static observational dataset rather than on discovering laws for dynamical systems, and 4) rather than a stochastic search over the space of functional forms, our proposed system leverages sophisticated techniques for enumerating expressions of increasing complexity (Larson and Van Cleemput, 2016 ). In our system, the human remains "in the loop" to evaluate the plausibility of suggested bounds.
Langley (2019) provides a review of past efforts in computational scientific discovery. Several frameworks have their origins in analyzing mass spectroscopy and other electrochemical data. Bacon (Langely et al, 1987 ) is a general framework for scientific discovery based on suggesting and executing a series of designed experiments. Tallorin et al (2018) proposed a method called POOL that uses Bayesian optimization and machine learning in an iterative fashion for experiments to discover peptide substrates for enzymes. Bacon and POOL both make recommendations regarding additional data to collect while our system assumes that a fixed dataset is provided that may or may not be the result of a designed experiment.
Interpretable AI. Precise definitions of "explainability" and "interpretability" are still being developed (Vilone and Longo, 2020; Lu et al, 2019; Fürnkranz et al, 2020) . According to the convention of Rudin (2019) , explainability is concerned with post-hoc analyses of black box models to create simple explanations of model behavior. Motivated by observed accuracies of deep learning models, work in this area includes identifying important features for prediction, building simple local models, conducting sensitivity analyses, and deriving prototype examples (Samek and Müller, 2019; Elton, 2020) . Rudin (2019) advocates instead the development of interpretable models where the mechanisms for predictions are simple relationships that are readily apparent to the investigator. Most recent work in this area is in the development of tree-based or rule-based models (e.g., (Wang et al, 2017; Rudin and Ertekin, 2018) ). Such models can be viewed as a series of nested if-then clauses. Different from these works, our proposed framework automates the discovery of nonlinear features. In addition, we provide a method for automatically discovering boolean relationships which is an alternative approach to rule-based models.
Automated feature generation. An important distinction of our work from prior work on automated feature generation is that prior work is focused on improving accuracy, sometimes at the expense of understandable features, and not on scientific discovery. Khurana et al (2018) propose a system that leverages reinforcement learning to search expression trees for predictive features. ExploreKit (Katz et al, 2016 ) is a framework for automatic feature engineering that combines features using basic arithmetic operations and then uses machine learning to rank and select those with high predictive ability. The Data Science Machine (Kanter and Veeramachaneni, 2015) automatically generates features for entities in relational databases with possible dependencies between tables followed by singular value decomposition. In none of these works is model transparency evaluated but rather only model performance.
Statistical methods. Traditional statistical methods for empirical model building (e.g. regression analysis) tend to focus on first-and second-order polynomial models; interaction terms up to a certain degree are often included. Em-pirical models are intended to provide adequate prediction performance while also providing a simple assessment of feature importance via model coefficients. Techniques such as all-subsets, stepwise selection, and regularization methods are commonly used to perform feature selection over model spaces of increasing complexity. However, domain knowledge is typically required for non-polynomial relationships (including reciprocal relationships). While nonlinear regression techniques exist (Seber and Wild, 2003; Song et al, 2014) , searching for such relationships is a departure from common practice in statistics. Our proposed framework provides a search over a much broader class of nonlinear functions.
The remainder of this paper proceeds as follows. In section 2, we illustrate the inability of traditional machine learning techniques to reveal nonlinear features to the user. Two illustrative examples are provided. Section 3 introduces our computational conjecturing framework for scientific discovery and automated feature engineering, and includes results for the examples from section 2. Section 4 concludes the article with a summary and discussion.
In this section, we provide two demonstrations where a "typical" knowledge discovery workflow fails to reveal important relationships among features. The workflow that we employ begins with simple statistical models and is followed by more complex models and modern machine learning methods.
Research on machine learning does, of course, lead to conjectured relationships between variables which are in turn used to make predictions of one or more variables in terms of others. A trained neural net, for instance, can be viewed as a black box representing a function which produces an output for every input in its domain. These functions are complex and of a different character than classical scientific laws: in particular, there is little hope of deriving these functions or relationships from simpler existing laws.
The first example is a simple case where a binary target variable is completely determined by the product of two continuous features in the dataset; i.e., the second-order interaction term completely defines the relationship.
Consider a dataset on residential real estate properties for sale obtained from https:\www.redfin.com. The machine learning goal is to predict whether a home with given attributes has a list price above or below $300,000.
This dataset includes both the price per square foot and total square footage along with eight additional features such as the number of bathrooms and bedrooms. The target (above vs. below) can obviously be determined by simply multiplying the price per square foot by square footage and setting a threshold. Thus, the interaction of price per square foot and square footage, hereafter called the active interaction, completely describes the relationship between the predictors and response. We investigate the ability of various statistical and machine learning techniques to leverage and reveal this relationship.
We begin with a typical statistical approach for understanding the relationship between the predictors and the target. A logistic regression model with only main effects produces an area under the curve (AUC) of over 0.98. Since the active interaction term is not included in the model but the model has high accuracy, an investigator might not even consider searching for interaction effects and therefore miss the fact that the product of two of the features determines the response. For the fitted model, the terms comprising the active interaction (i.e., price per square foot and square footage) are statistically significant predictors (using a standard significance level of α = 0.05), along with several other variables.
The fitting algorithm for logistic regression fails to converge if all 10 2 = 45 two-way interaction terms are added to the main effects model. If predictors are first screened based on the fitted model with only main effects, the algorithm for the model with the seven statistically significant predictors and their 21 interaction terms fails to converge as well.
When LASSO (Friedman et al, 2010 ) is applied to centered and scaled data for the model with all main effects and 45 interaction terms, the active interaction term shrinks to zero with nearly two-thirds of the solution path remaining. Further, the LASSO model that minimizes cross-validated error does not contain the active interaction term. Now consider a non-standard ensemble method based on generalized linear models that is designed to leverage interaction terms up to a specified degree called randomGLM (Song et al, 2014; Song and Langfelder, 2013) . The final model is based on a bootstrap aggregation across random samples of observations and features (bags). For this data, the final model produces an AUC of over 0.99 with the active interaction term occurring in more than half of the bags. Over 15 other interaction terms occur in at least one bag as well.
Finally, we consider several machine learning methods including classification trees, support vector machines (SVM), neural networks (NN), and random forests (RF). These methods are designed to leverage, but not explicitly discover, complex relationships among predictors. NN and RF have easily-accessible methods for evaluating the contribution of individual variables to prediction accuracy, but do not reveal nonlinear relationships among variables.
Models were fit in R using the libraries rpart ( (Liaw and Wiener, 2002) (RF). Tuning was performed using caret (Kuhn, 2020) . The code and results are available on Github https://github.com/nvcleemp/conjecturing. The classification tree model returns high importance scores for the variables in the active interaction along with several other variables and the AUC is 0.88. SVM with a linear kernel achieves an AUC of over 0.98, mirroring the results observed for logistic regression with only main effects. SVM with a degree two polynomial kernel can leverage interaction terms for prediction; the resulting model also has an AUC of over 0.93. SVM with a radial basis function kernel achieves an AUC of 0.89.
An NN model built using default settings with tuning for the number of hidden units and the weight decay achieved an AUC of over 0.99. Price per square foot and square footage are identified as the most important predictors by Garson's algorithm (Garson, 1991; Beck, 2018) . Other variables had positive relative importance values as well. An RF model achieves the highest AUC among the models we tested, above 0.99. As with the classification tree model, price per square foot and square footage have the highest variable importance scores among all predictors, but it is not possible to determine from the algorithm that they should be multiplied. Figure 1 contains receiver operating characteristic (ROC) curves for the various methods tested. The ROC for Random GLM is the best, followed by those for NN and RF. Logistic regression, LASSO, and SVM with a linear kernel had slightly worse ROCs. The implementations of SVM with the polynomial kernel and SVM with the radial basis function kernel failed to provide class probabilities for test instances, so the AUC is only based on one false positive rate.
From this example, it is clear that some of the machine learning methods are able to achieve high rates of prediction accuracy and some can identify the terms of the active interaction term as important, but none are helpful to the investigator for discovering that the terms should be multiplied.
In this example, a continuous target variable is determined by a more complex nonlinear relationship with three continuous predictors. We again suppose that the investigator is unaware of the nonlinear relationship and applies a standard analysis using statistical and machine learning tools. Although some of the machine learning methods are able to leverage the patterns to generate good predictions, the analysis fails to indicate the true nonlinear relationship.
Consider measurements including the masses of two objects m 1 and m 2 , their distance r, and the gravitational force between them F . The goal is to recover the dependence of F on m 1 , m 2 , and r, or
where k is the gravitational constant. Following the demonstration by Langely et al (1987) , we create a fictional dataset using a predefined value for k that is a random number between 0 and 1. For our illustrative example, we generated 1,000 points with k = 0.057098. As in the real estate valuation example, the target F can be determined by evaluating a nonlinear expression of the input variables m 1 , m 2 , and r. We again investigate the ability of various statistical and machine learning techniques to leverage and reveal this relationship. The ability of each method to predict the force F for a holdout dataset is also evaluated.
A regression model with only main effects produces a model where each term is identified as statistically significant at a significance level of 0.05. We obtain an R 2 of 0.090 indicating a weak linear relationship. Although it is helpful that the regression analysis informs the user that the relationship is more complex than specified, there is nothing to suggest the true functional form.
A second-order regression model identifies as significant the terms r 2 , m 1 m 2 , m 1 r, and m 2 r. The first term is inversely proportional to F and the second term is directly proportional to F , but the true functional relationship still remains undiscovered.
The RF model performed best with the lowest prediction error. The three terms m 1 , m 2 , and r were identified as important via variable importance measures, but the nonlinear relationship is hidden to the investigator. The regression tree, SVM, and NN models also outperformed the regression models in terms of prediction but by their design, do not reveal the relationship. Figure 2 presents the root mean square error (RMSE) for test observations for the various methods tested. Only RF performed noticeably better than the linear model with only main effects with an RMSE less than half of the other methods. We suspect that the reason for this disparity is the fact that the data are not normalized.
As with the real estate example, a machine learning method can leverage the nonlinear relationship in the data but cannot present the relationship to the investigator. In the next section, we propose a framework for discovering the nonlinear relationships in both cases.
We now describe a framework that leverages a conjecturing algorithm to discover nonlinear and boolean feature relationships that can enhance understanding and predictive ability. We assume that data x i ∈ R d , i = 1, . . . , n are given. Each of the n points correspond to examples for the conjecturing algorithm and the d features are invariants. There are two versions of the conjecturing algorithm (Conjecturing): one that produces invariant relations (Conjecturing-INV) and one that produces property relations (Conjecturing-PROP).
An invariant is a numerical attribute of an example. Conjecturing-PROP searches for bounds (upper or lower) on an invariant of interest that are functions of the other invariants. The bounds must hold for all examples. An invariant relation is retained only if the upper/lower bound is tighter than all previouslyfound bounds for at least one example. Figure 4 displays (a) upper bounds and (b) lower bounds derived for test instances for the gravity data. The gray curves correspond to bounds, and each must be the best on at least one training example in order to be retained.
A property is a boolean attribute of an example. The property relations version searches for logical expressions that are true for every example. One of the properties is specified as the property of interest. The conjectures that are produced are conditional statements that are sufficient conditions for the property of interest complexity of one. Expressions are formed from the invariants/properties by applying operators and are represented using an expression tree ( Figure 5 ). The complexity of an expression is the number of nodes in the expression tree.
Our implementation of the conjecturing algorithm on a computer with an Intel Core i5-4210U 1.7 GHz processor and 4 GB RAM generates over 7.2 billion expressions using 20 invariants/properties in one minute. At one minute, expressions of up to complexity six are generated. Figure 5 (a) depicts the expression tree for an expression of complexity five involving features from the real estate data. Figure 5 (b) depicts an expression tree for an expression of complexity eight for the formula for gravity, assuming that the gravitational constant k is known. By default, the conjecturing algorithm includes operators such as adding or subtracting 1, dividing or multiplying by 2, and powers of 10. A gravitational constant of 0.057 could be represented as (+1 + 1 + 1 + 1 + 1)/10 −1−1 + (+1 + 1 + 1 + 1 + 1)/10 −1−1−1 which by itself has complexity 19. We will use the conjecturing framework to recover the functional form of gravity without k which is an expression of complexity 7.
More details on Conjecturing are included in the Appendix. For the computational experiments described in the remainder of the paper, we use an implementation of Conjecturing that is open-source and operates using SageMath (https://sagemath.org). The program, examples, and set-up instructions are available at: http://nvcleemp.github.io/conjecturing/.
Our proposed framework leverages the invariant version (Conjecturing-INV) and then the property version (Conjecturing-PROP) of the conjecturing algorithm (Algorithm 1). For a classification problem, the framework is designed to produce new binary features that indicate whether nonlinear patterns are satisfied by the original feature values for each observation. These new features can augment the original features in a classification model, capturing nonlinear patterns while maintaining interpretability.
We assume that a dataset with class labels (x i , y i ) ∈ R m × Y, i = 1, . . . , n is given. We also assume that categorical features including the class y have been encoded as a series of binary features. For each class y ∈ Y, the algorithm discovers inequalities that are satisfied by each observation in the class (Steps 4-14) . These inequalities are converted to properties of the form "if the inequality is satisfied, then True; False, otherwise" (Step 12). These new properties are combined with the original binary features in the data (Step 13). The properties from across all classes are pooled together and the observations belonging to all classes are pooled together as examples and then, for each class y, the property version of Conjecturing is applied to discover sufficient conditions of the form "if the logical expression is satisfied, then the observation belongs to the class" (Steps 15-19) . These new properties are then returned and can be incorporated in a classification model (e.g. a logistic regression model).
We now provide further details on Algorithm 1 using the real estate valuation case from Section 2.1 as an illustrative example. First, we binarize the categorical attribute propertyType into binary features condo, mobileHome, singleFamily, townhouse, multiFamily2-4Unit, multifFamily5PlusUnit, and Other. We also add an attribute that is a constant value of 300,000 for each observation because it is the price cutoff that determines the target and call it threeHundredK. 
Also included are less-interpretable bounds such as:
hoaPerMonth ≤ 10 2×bathrooms + squareFootage
There were also several bounds discovered that are close approximations of the relationship present in the active interaction term, including squareFootage ≤ threeHundredK/pricePerSquareFoot + bathrooms (7) squareFootage ≤ threeHundredK/pricePerSquareFoot + bedrooms (8) squareFootage ≤ threeHundredK/pricePerSquareFoot + daysOnMarket (9) squareFootage ≤ threeHundredK/(pricePerSquareFoot − 1) − 1 (10)
pricePerSquareFoot ≤ −threeHundredK/(bedrooms − squareFootage)
pricePerSquareFoot ≤ threeHundredK/squareFootage (12) threeHundredK ≥ −(bathrooms − squareFootage) × pricePerSquareFoot (13) For houses with target value above, there are 1,457 bounds derived including a mix of simple relations and less intuitive relations. Also included are the following three relations that are nearly identical to the active interaction relation:
pricePerSquareFoot ≥ threeHundredK/squareFootage + 1 (15) threeHundredK ≤ (pricePerSquareFoot + 1) × squareFootage
The resulting invariant relations are pooled together (Step 10). The invariant relations are encoded as properties (Step 12). The original binary features from the data are also encoded as properties for a total of 1, 280 + 1, 457 + 7 = 2, 744 properties. Examples of encoded properties from the invariant relations are:
(squareFootage ? ≥ threeHundredK/(pricePerSquareFoot + 1)).
The properties generated across the classes {below, above} are collected in a set Π along with the class label and the seven original binary features (Step 13).
For each class, apply the property version (Conjecturing-PROP) of Conjecturing to the properties Π with the entire training data observations serving as the examples E and class membership as the target property (Steps 15-19) . The result is a set of properties that are sufficient conditions for class membership.
Conjecturing-PROP returns only two properties. They both approximate the underlying active interaction.
An inspection of the data reveals that for some of the houses, there is some rounding error when comparing the price to the square footage multiplied by the price per square foot. The conjecturing algorithm compensates by using invariants as error terms. In the first property, the error term is bathrooms×pricePerSquareFoot. In the second property, the error term is squareFootage + 1.
When these properties are applied as classification rules, they produce no error on the training data. The first property misclassifies 37 of 30,156 houses in the test data for an accuracy of 0.999. The second property misclassifies 26 houses. The misclassified houses are due to rounding error and miscoding of data. For example, one house in the test data is listed as having 31,248 bathrooms and another is listed as having a price of $459.
Despite the noise and rounding error in the data, the conjecturing framework was able to recover the active interaction term and can be used as classifiers with near-perfect accuracy.
Our proposed framework for a continuous target variable is simpler than for classification. Only the invariant version (Conjecturing-INV) of the conjecturing method is used.
For a regression problem, the framework produces nonlinear functions of the original variables that can be used to predict the response. The functions can then be used as features in a regression model.
We assume that a dataset (x i , y i ) ∈ R m × R is given, where the y i are values for the response variable. Invariant conjectures are generated that provide upper and lower bounds on the response. These conjectures are the nonlinear functions that can be used as new features.
Consider the gravity case from Section 2.2. First, we center the values for force F so that the intercept for any fitted model is zero (Step 1). The invariants are F , m 1 , m 2 , and r (Step 3). The examples are the observations in the data (Step 4).
The conjecturing framework is not designed to recover constants such as the gravitational constant k. In general, for a functional relationship with a constant k such that 0 < k < 1, the expression without the constant provides a lower bound for the response and the reciprocal expression provides an upper bound. In cases where the constant is not between 0 and 1, the converse is true. For this example, F = km 1 m 2 /r 2 with 0 < k < 1 and so m 1 m 2 /r 2 ≤ F ≤ r 2 /m 1 m 2 . The conjecturing framework can potentially recover the bounds F ≥ m 1 m 2 /r 2 or F ≤ r 2 /m 1 m 2 .
For our example, Conjecturing-INV returns 19 upper bounds (Step 5) and 24 lower bounds (Step 6) for F . Among the upper bounds is F ≤ m 1 m 2 /r 2 , which approximates the true gravity relationship used to generate the data. The bound does not include the constant k. Other bounds generated by Conjecturing-Algorithm 2 Conjecturing framework for regression models. Input: Data (x i , y i ) ∈ R m × R, i = 1, . . . , n. Output: A set of invariant relations R. 
Eight of the upper bounds and 15 of lower bounds for F are depicted in Figure 3 . The upper bound m 1 m 2 /r 2 in Figure 3 (a) is blue, while the true value km 1 m 2 /r 2 is gold. In this example, the proposed conjecturing framework recovers the true nonlinear relationship while regression models and machine learning methods fail to do so.
In this section, we demonstrate the proposed framework on synthetic patient-level COVID-19 data that was provided as part of the VHA Innovation Ecosystem and precisionFDA COVID-19 Risk Factor Modeling Challenge (https://precision. fda.gov/challenges/11/view). The data include synthetic veteran patient health records including medical encounters, conditions, medications, and procedures.
The goal of the challenge was to predict COVID-19 outcomes including alive/deceased status and whether a patient would require ICU care or not. Predictions were based on information obtained through December 31, 2019. In the training data, we drop all information pertaining to events on or after January 1, 2020 and drop subjects who died before January 1, 2020.
We generate a number of features for each patient based on the data, described in the Appendix. In addition, for each reported allergy, device, immunization, procedure, and discretely-measured observation we create an indicator variable to serve as a property. In total, we use 309 invariants and 362 properties. We generate conjectures for three randomly-selected sets of 10 subjects for upper and lower bounds for each invariant, and for each outcome (alive/deceased, ICU/no-ICU). These bounds are used as properties for Conjecturing-PROP with samples of 10 subjects for each outcome, and the process is repeated for 100 iterations. Conjectures are generated for both outcomes.
Alive/Deceased Status. Our analysis produces six conjectures as sufficient conditions for alive status and five conditions for deceased status as shown in Table  1 . 
Among those with COVID-19 in the data, 5,568 (7.6%) have a status of deceased, and 68,129 (92.4%) are alive. There are four single-factor conditions: having ever had a pneumococcal polysaccharide vaccine, alcoholism, hyperlipidemia, and osteoprosis. For patients with these conditions, the percentage of those dead from COVID-19 is larger than the 7.6% in the general population (8.5%, 9.3%, 12.1%, 16.4%, respectively).
Among those with the last sufficient condition,
16.1% died from COVID-19, so that subjects satisfying this condition died at a rate more than twice as high as for the general population. This condition involves the lifetime healthcare expenses for a subject along with the longitude of their home address. Figure 6 (a) contains a plot of subjects by location and COVID-19 death status and Figure 6 (b) contains a plot of subjects by location and lifetime healthcare expenses. It is immediately apparent that all subjects in the dataset reside in Massachusetts. Approximately the eastern third of Massachusetts is comprised of Boston and its metropolitan area. The population density, along with higher numbers of COVID-19 deaths, are notable around the city center. The population density of Springfield in the southwest is also apparent.
The bound indicates that, for a given level of healthcare expenses, those residing further to the east are more likely to have a status of deceased. Indeed, most of the COVID-19 deaths are located in the east of the state. The property seems to suggest that risk factors include living in the densely-populated east and living in the west with higher-than-typical lifetime healthcare expenses. The bound itself provides an easily-calculable quantity that can help to identify high-risk patients. Five of six of the sufficient conditions for alive status are true for more than 92.4% of the alive COVID-19 patients.
ICU Status. The framework produces seven conjectures as conditions for not requiring ICU care and six conjectures as conditions for requiring ICU care among patients with COVID-19 (see Table 2 ). Among those with COVID-19, 4,981 (7.2%) require ICU care while 68,716 do not (92.8%).
There are three single-factor conditions: hyperglycemia disorder, coronary heart disease, and former smoker status. Hyperglycemia is a symptom of diabetes which, along with coronary heart disease and smoking, is recognized as a risk factor for severe illness from COVID-19 (https://www.cdc.gov/coronavirus/2019-ncov/ need-extra-precautions/people-with-medical-conditions.html, accessed October 19th, 2020) . Patients with these conditions require ICU care at higher rates than the general population of COVID-19 patients (9.8%, 12.9%, and 8.5%, respectively). The conjectures involving more attributes are not true for more than 7.2% of the population and so are not as accurate as choosing a patient at random.
Five of the seven sufficient conditions for non-ICU care are true for more than 92.8% of COVID-19 patients that did not receive ICU care. Among these, 98.3% of those satisfying the left-hand side of the conjecture healthcareCoverage > (encountersLifetimeTotalCost − 1) × age → no-ICU, do not require ICU care. This conjecture suggests an interaction between the cost of medical care, age, and healthcare coverage. 
We have demonstrated that automated search for conjectured feature-relations can enhance existing learning algorithms. The discovery of these kinds of feature relationships can also initiate new collaboration with domain scientists and lead to new scientific knowledge. Our proposed Conjecturing framework was able to recover a hidden interaction between price per square foot, square footage, and price in real estate data that leads to improved classification performance. The framework also recovered the functional form for gravity with only the measured force, masses, and distance. Using recent patient-level COVID-19 data, the framework produced conjectures that provide insight into the risk of a need for ICU care and the risk of death.
The current version of Conjecturing requires that conjectures are true for every example. Future research will adjust the algorithm to better handle noisy data by generating conjectures that do not necessarily hold for all examples.
If the Conjecturing framework can provide the functional relationship without the constant of proportionality, the constant can be determined using regression with the original data. Suppose that the conjecturing framework indicates a relationship between the response y and predictors x of the form y ≤ b 1 f (x) for an unknown constant b 1 . A regression model can be fit of the formŷ = b 0 + b 1 f (x) using the data (x i , y i ), i = 1, . . . , n.
Another area for potential research involves the so-called p >> n problem. That is, if the number of features is larger than the number of observations, then there are insufficient degrees of freedom to estimate a linear model with all p features or any more complex model. In such situations, feature and/or model selection tools are needed to search over potentially large model spaces. As desired model complexity increases (e.g. consideration of interaction terms), searching over such large model spaces can become computationally prohibitive. For instance, suppose an investigator seeks to identify a model by selecting the "best" subset from among 10 features and their associated 45 two-way interactions. In this example, simply considering models with only 10 variables requires searching over a model space larger than 2.9 billion. Future research will investigate the ability of the Conjecturing framework to simplify model spaces and hence, provide a mechanism for a more expeditious search of plausible models.
The algorithm we have used to conjecture feature relations is an adaptation of an algorithm that was originally designed to conjecture invariant relations for mathematical objects.
In the case of graph theory, for instance, the independence number of a graph, the largest number of vertices of the graph for which no pair of vertices is contained in an edge, is a widely-studied NP-hard-to-compute graph invariant. Bounds for this invariant are of interest for both theoretical and practical reasons: bounds for instance can be used to help speed up computation of this number. The program we use was originally designed, for instance, to compute upper and lower bounds for the independence number of a graph.
The program is given the following inputs: an invariant of interest, a bound type (upper or lower) some number of example graphs, and some number of other graph invariants. The program must be able to compute the value of each invariant for each given graph (or at least be provided with their values). Furthermore, the program will produce expressions involving some number of provided unary and binary operators: these are built-in to the program, and not ordinarily modified by the user. Currently these include unary operations like "+1", squaring, squarerooting, and division by 2, and binary operations including "+", "×", and "−". These can be changed by the user, if desired, in the code.
The program then generates expressions, of progressively greater complexity. The simplest expressions, of complexity 1, are the invariants themselves. The next simplest are unary operators applies to the given invariants; these have complexity 2. In general, the complexity of an expression beginning with a unary operator is one plus the complexity of the remaining expression, while the complexity of an expression beginning with a binary operator is one plus the complexities of the two sub-expressions. In this way all expressions, with progressively greater expressions, can be generated. These expressions are represented by a tree in the program: complexity is defined in terms of the order of the representing tree. Exhaustive generation of all possible expressions of a given complexity corresponds to a tree generation problem.
The program then produces expressions of the form "(invariant of interest) (bound-type) (expression generated from other invariants)". These become statements, true or false, when they are interpreted as being quantified over some domain. The program will continue to generate more and more statements until it reaches a program-defined stopping point: this stopping point may be, variously, when the program can no longer make a better conjecture (in a well-defined sense) or when some timing condition has been reached.
Each produced expression is tested for truth with respect to the examples known to the program. This is the truth test. Secondly, each produced statement must also be non-dominant. A statement is non-dominant if there is at least one object known to the program where the statement gives a better bound for the object than any of the previously accepted statements. For instance, if the program is generating upper bounds for the independence number α of a graph, and knows graphs G 1 , . . . , G k and is considering a potential conjecture C of the form α ≤ β; if C is true for all the graphs known to the program then, and if there is some graph G with β(G) less than the bounds given by all the other previously accepted conjectures then C is non-dominant and is accepted and added to the store of conjectures.
This non-dominance test is Siemion Fajtlowicz's "Dalmatian heuristic" (Fajtlowicz, 1995) and has been used to produce a large number of conjectures that have been investigated by mathematical researchers and have led to publications.
We now describe the Dalmatian heuristic more formally. Let O 1 , . . . , O n be examples of objects of a given type. Let α 1 , . . . , α k be real number invariants. And let α be an invariant for which conjectured upper or lower bounds are of interest. We form a stream of algebraic functions of the invariants: α 1 + α 2 , √ α 1 , α 1 α 3 , (α 2 + α 4 ) 2 , etc. These expressions can then be used to form conjectured bounds for α. If we are interested in upper bounds for α, we can form the inequalities α ≤ α 1 + α 2 , α ≤ √ α 1 , α ≤ α 1 α 3 , α ≤ (α 2 + α 4 ) 2 , etc. These inequalities can be interpreted as being true for all the objects of the given type. The inequality α ≤ α 1 + α 2 can be interpreted, for instance, as, "For every object O, α(O) ≤ α 1 (O) + α 2 (O)." A conjectured upper bound u is only added to the database of conjectures if the bound passes the following two tests. That is, the candidate conjecture would give a better bound for α(O) than any previously conjectured (upper) bound.
The general approach to generating conjectures is as follows.
1. Produce a stream of inequalities with evaluable functions of the invariants on each side of the inequality symbol. Some of these will pass the truth and nondominance tests and be stored as conjectures. 2. Initialize an initial collection of objects. These can be as few as one. 3. Generate conjectures that are true for all stored objects and non-dominant with respect to these objects and the previously stored conjectures. Pass each generated statement to the truth and non-dominance tests. The program needs a stopping condition. The best case is that, for each object, there is at least one conjecture that gives the exact value for the object. In this case there is no possibility of improving the current conjectures-in the sense that no other conjectures can make better predictions about the values of the existing objects-exact predictive power for all objects has been achieved. In the case where this natural stopping condition is never attained, another stopping condition will be required. One possibility is to simply stop making conjectures after some hardcoded or user-specified time. 4. Remove dominated conjectures. After a conjecture is added to the store of conjectures, it may be the case that another conjecture in the store is no longer non-dominant. If conjectured upper bounds (for example) for an invariant α are being generated then a conjectured bound α i in the conjectures store is non-dominant, with respect to the stored objects, if and only if there is an object O such that α i (O) < min{α j (O) : j = i}, that is, if and only if, there is an object O where the bound gives a better predicted value for α(O) than any other conjectured bound does. Dominated conjectures are then removed.
If dominated conjectures are removed whenever non-dominant conjectures are added to the conjectures store, the number of conjectures (of any particular form) cannot exceed the number of example objects. This puts an important limit on the number of produced conjectures.
Feature Name Definition healthcareExpenses
The total lifetime cost of healthcare to the patient healthcareCoverage
The total lifetime cost of healthcare services that were covered by payers latitude
Latitude of patient's home address longitude
Longitude of patient's home address age Current age of patient 