Optimal design of experiments by combining coarse and fine measurements

TL;DR

This paper introduces a strategy combining coarse and fine measurements to efficiently build accurate predictive models, reducing costly high-resolution data needs by leveraging abundant categorical data.

Contribution

It presents a novel approach inspired by statistical physics that integrates coarse and fine measurements for improved model accuracy.

Findings

Effective combination of coarse and fine data improves prediction accuracy.

Method applied successfully to predict molecular properties.

Reduces need for extensive high-resolution measurements.

Abstract

In many contexts it is extremely costly to perform enough high quality experimental measurements to accurately parameterize a predictive quantitative model. However, it is often much easier to carry out large numbers of experiments that indicate whether each sample is above or below a given threshold. Can many such categorical or "coarse" measurements be combined with a much smaller number of high resolution or "fine" measurements to yield accurate models? Here, we demonstrate an intuitive strategy, inspired by statistical physics, wherein the coarse measurements are used to identify the salient features of the data, while the fine measurements determine the relative importance of these features. A linear model is inferred from the fine measurements, augmented by a quadratic term that captures the correlation structure of the coarse data. We illustrate our strategy by considering the problems of predicting the antimalarial potency and aqueous solubility of small organic molecules from their 2D molecular structure.

Equations

$$C_{\pm }$$

$$=\sum _{i=1}^{N_{\pm }}{\lambda }_i^{\pm }{\mathbf{u}}_i^{\pm }\otimes {\mathbf{u}}_i^{\pm },$$

$$y_i={\mathbf{h}}^T{\mathbf{f}}_i+{\mathbf{f}}_i^{T}J{\mathbf{f}}_i+{ϵ}_i.$$

$$J=\sum _{k=1}^{{\hat{p}}_{+}}{c}_k^{+}{\mathbf{u}}_k^{+}\otimes {\mathbf{u}}_k^{+}+\sum _{k=1}^{{\hat{p}}_{-}}{c}_k^{-}{\mathbf{u}}_k^{-}\otimes {\mathbf{u}}_k^{-}.$$

$$J=\sum _{i=1}^m{E}_i{\mathbf{\zeta }}_i\otimes {\mathbf{\zeta }}_i.$$

$$\rho (\lambda )=\frac{\sqrt{{\left[{\left(1+\sqrt{\gamma }\right)}^2-\lambda \right]}_{+}{\left[\lambda -{\left(1-\sqrt{\gamma }\right)}^2)\right]}_{+}}}{2\pi \gamma \lambda }$$

$$H(\mathbf{f})=\sum _i{E}_i\exp \left(-E_i^2(\mathbf{f}\cdot {\mathbf{\zeta }}_i{)}^2\right).$$

Full text

A large class of scientific questions asks whether dependent variables can be accurately predicted by using training data to learn the parameters of quantitative models. Classical statistics shows that this is possible if sufficiently many high resolution measurements are available, though the cost of performing these experiments can be prohibitive. On the other hand, in many settings, it can be straightforward to evaluate whether a measurement is above or below a certain threshold, raising the question of how such measurements can be incorporated into the modelling framework.

Examples abound in disparate fields. For instance, predicting the solubility of organic molecules is a fundamental challenge in physical chemistry Llinas et al. (2008). Although accurate measurements are extremely difficult to obtain Palmer and Mitchell (2014), determining whether a molecule is soluble at a particular concentration is comparatively simple. Similarly, in drug discovery, biochemical assays that determine whether a molecule binds to a given receptor are much simpler than measuring protein-ligand binding affinity Malo et al. (2006). In protein biophysics, a key challenge is to predict the effect of amino acid changes on protein phenotype. Here, threshold measurements are naturally provided by homologous sequences from the same protein family Morcos et al. (2014); Figliuzzi et al. (2015); Hopf et al. (2017); Barrat-Charlaix et al. (2016); Levy et al. (2017). In contrast, experimentally measuring the phenotypic change is much more difficult. A related problem is to predict the viral fitness landscape given HIV sequences obtained from patients; again collecting patient samples is much easier than measuring fitness directly Shekhar et al. (2013); Ferguson et al. (2013). In single-cell RNA sequencing, decomposition methods that extract the correlation structure of shallow gene expression measurements is an ongoing challenge Dixit et al. (2016); Buettner et al. (2016).

Despite the ubiquity of this problem, to our knowledge there is no principled method for combining numerous binary/categorical (“coarse”) measurements with fewer quantitative (“fine”) measurements to produce a predictive model. Although regression approaches can account for a prior estimate of sample error Bishop (2007), this is not the same as combining two qualitatively distinct forms of data to build a more accurate model.

In this Letter, we introduce an intuitive method that combines coarse and fine measurements. The coarse measurements provide sets of labelled samples – those data above and below some threshold – and the proposed method extracts features from the correlations of the variables in each set. Their contribution to the dependent variable is then determined by using the fine measurements to build a regression model for these features. Our model augments a quantitative linear model with a quadratic term which captures the correlation structure extracted from the coarse data. We illustrate our approach by applying it to solubility prediction, and interpret the approach in the context of the Ising model.

To fix ideas, we assume each sample is characterized by a vector of $p$ properties $\mathbf{f}_{i}\in\mathbb{R}^{p}$. The binary data indicates that $N_{+}$ ($N_{-}$) samples are above (below) some threshold. In addition, we are given $\mathbf{y}\in\mathbb{R}^{M}$, quantitative measurements for $M$ additional samples. These measurements could be binding affinity, solubility etc. We construct matrices $R_{\pm}\in\mathbb{R}^{N_{\pm}\times p}$ for samples above/below the threshold, with columns of $R_{+}$ and $R_{-}$ normalized separately to have zero mean and unit variance. Intuitively, if there are combinations of the $p$ properties that are always present in either sample set, then these properties should be good predictors of the measurement. Such persistent correlations can be identified from the eigendecomposition of each sample covariance matrix $C_{\pm}$

[TABLE]

where $\{\lambda^{\pm}_{i}\}$, $\{\mathbf{u}^{\pm}_{i}\}$ are the eigenvalues and eigenvectors (note we perform separate eigendecompositions for the two matrices $C_{\pm}$). Each $\mathbf{u}^{\pm}_{i}$ identifies a particular combination of the $p$ properties, explaining a fraction $\lambda^{\pm}_{i}/\sum_{i}\lambda^{\pm}_{i}$ of the variance Bishop (2007). Each matrix $C_{\pm}$ is an unbiased estimator of the corresponding true covariance matrix. The quality of this estimator depends on data sampling. For example, one may inadvertently assay certain samples (easy to obtain, measure etc.), which could distort the estimator by causing an eigenvector with large eigenvalue to be localized on features common to these samples, even though they do not predict the output variable. For protein sequences, a natural source of such spurious correlations is phylogeny Dutheil (2012); Obermayer and Levine (2014).

Here we propose that whereas the eigenvectors $\mathbf{u}^{\pm}_{i}$ reliably identify data features, their significance as estimated by the corresponding eigenvalues $\lambda^{\pm}_{i}$ can be severely corrupted by imperfect sampling. Later we justify this ansatz with ideas from statistical physics, and show that this characterization applies to a large class of problems. This ansatz suggests a strategy to mitigate the corruption by using the additional quantitative measurements to determine the significance of each feature. We posit a general quadratic model

[TABLE]

Here $\mathbf{h}$ is the variable-specific effect, $J$ captures the coupling between variables, and $\epsilon_{i}\sim N(0,\sigma)$ models random error. There are $p$ parameters in $\mathbf{h}$ and $p(p-1)/2$ parameters in $J$. If one had $M\gg p(p+1)/2$ quantitative measurements, these parameters could be estimated using linear least squares regression. However, it is costly to perform many detailed measurements, so we turn instead to the matrices $C_{\pm}$. We pose the ansatz

[TABLE]

Here $\hat{p}_{\pm}\leq p$, since some eigenvectors will reflect noise due to finite sampling Laloux et al. (1999); Plerou et al. (1999); Bai and Silverstein (2010). Our ansatz reflects the hypothesis that the eigendecomposition of $C_{\pm}$ captures variable-variable correlations. If the number of samples is much smaller than the number of variables, random matrix theory provides a rigorous way to determine $\hat{p}_{\pm}$ Laloux et al. (1999); Plerou et al. (1999); Bai and Silverstein (2010); Bun et al. (2016); Lee et al. (2016); Bun et al. (2017); this case will be discussed in detail later. Relaxing this assumption, we include all eigenvectors and determine the parameters $\mathbf{h}$, ${c}^{+}$ and ${c}^{-}$ by regressing against the the few quantitative measurements available. We note that the ansatz (3) reduces the number of variables to $p+\hat{p}_{+}+\hat{p}_{-}$. In the case where the coarse measurements yield multiple categories (or a single category), our method generalizes by forming separate correlation matrices for each category, and positing that $J$ is a sum of the outer product of all eigenvectors with coefficients determined by regressing against the quantitative data. Our method generalises to Generalised Linear Models with a link function on the right hand side of Equation (2).

To illustrate our approach, we consider two examples: predicting the potency of chemicals against malaria and the equilibrium aqueous solubility of molecules.

Antimalarials – Developing accurate models that can rank the potency of a library of compounds against a target is an important unsolved challenge in drug discovery. We consider a published antimalarial screening campaign Guiguemde et al. (2010): binary but high throughput assays reported 1528 active compounds against malaria, lower throughput but quantitative assays measured the potency ($\mathrm{pIC}_{50}$) of only 1189 compounds Guiguemde et al. (2010); Riniker et al. (2017). Figure 1A shows that by combining binary and quantitative measurements, an hitherto unattempted strategy, an order of magnitude less quantitative measurements could have been performed to yield a model with similar predictive accuracy (c.f. Supplementary Information showing the same result for the Pearson correlation coefficient). Moreover, the model with coarse measurements clearly outperforms the linear model without the coarse measurements and the “null” quadratic model where the vectors $\mathbf{u}^{\pm}_{k}$ are random orthogonal vectors (i.e. random $R_{\pm}$). In the Supplemental Information we show that our model also outperforms direct quadratic regression. The compounds are described using the 1024-bit Morgan6 Fingerprint Rogers and Hahn (2010) generated with rdKit rdk .

Solubility – Predicting the aqueous solubility of molecules is a fundamental problem in physical chemistry important to a plethora of chemical industries. However, accurate solubility assays are low throughput ($\sim$ 1 hour/compound Box et al. (2006)). Figure 1B shows that one could obtain an accurate solubility model ($r^{2}=0.85$, MAE = 0.61) if one were to combine the outcome of a coarse solubility assay that could only tell whether a compound is soluble ($<10^{-4}~{}\mathrm{mol/L}$) or not ($>10^{-2}~{}\mathrm{mol/L}$) with much fewer quantitive solubility data. We use a standard dataset of the solubility of 1144 organic molecules Delaney (2004), and describe the molecule by concatenating the Avalon Fingerprint Gedeck et al. (2006), the MACCS Fingerprint Durant et al. (2002), and the 1024-bit Morgan6 Fingerprint Rogers and Hahn (2010). Our result compares favourably with other models that also use binary molecular fingerprints, e.g. kernel partial least squares regression achieves $r^{2}=0.83$ Zhou et al. (2008).

To understand why our heuristic strategy is successful, we consider a model problem where data is generated according to Eqn. (2), which is the maximum entropy model Schneidman et al. (2006); Cocco et al. (2009); Lee et al. (2015), analogous to the Ising model. We thus interpret the dependent variable as an “energy”, noting that the logarithm of the solubility is proportional to the solvation energy. The interaction matrix $J$ can be decomposed into a sum of outer products of eigenvectors $\boldsymbol{\zeta}_{i}$ (Hopfield patterns Hopfield (1982)), and eigenvalues $E_{i}$ (Hopfield energies) as

[TABLE]

Furthermore, to model the binary features used in solubility prediction, we make the assumption that the independent variable is a vector of $\pm 1$.

To simulate binary measurements we randomly draw samples from the uniform distribution, evaluate Eqn. (2) to determine the energy of each sample, and retain those samples that fall below a certain energy. Consider an interaction matrix $J$ with $p=100$, and $m=3$ randomly generated patterns. To fix ideas, henceforth let all patterns be attractive with $(E_{1},E_{2},E_{3})=(-30,-25,-20)$, $\mathbf{h}=0$ and $\epsilon_{i}=0$. Using this model, we generate 5000 random vectors, consider the $N=500$ samples with lowest energy as above threshold samples, and compute the eigendecomposition of the resulting correlation matrix.

Figure 2A shows that the eigenvalue distribution of the sample correlation matrix $C_{+}$ follows the Marčenko-Pastur distribution expected for a random matrix Marčenko and Pastur (1967),

[TABLE]

where $(\cdot)_{+}=\mathrm{max}(\cdot,0)$, $\gamma=p/N$, with the exception of three distinct eigenvalues. Figure 2B shows that their corresponding eigenvectors are indeed the Hopfield patterns that we put in. Therefore, the large eigenvectors of $C_{+}$ correspond to eigenvectors of $J$. Note that the random matrix theory framework applies because $m\ll p$, i.e. the signal is low rank compared to the noise. If the signal was high rank, all eigenvectors should be included and their significance determined by regression against fine measurements, as in the examples discussed above.

Turning to eigenvalues, in this model, which features uniform sampling, we find that the eigenvalues are proportional to the Hopfield energy $E_{i}$. This allows us to “clean” the correlation matrix, by using the $q$ eigenvectors above the Marčenko-Pastur threshold to construct a rank $q$ approximation $J_{\mathrm{MP}}$ of the correlation matrix (here $q=3$). Figure 2C shows that $J_{\mathrm{MP}}$ accurately reconstructs $J$, and allows accurate prediction of the energy of particular states (Figure 2D). Analogously, the eigendecomposition of $C_{-}$ allows one to recover repulsive patterns with positive Hopfield energies (see Supplemental Information).

Since the Hopfield patterns are energy minima, taken together Figure 2A-D imply that the probability of the system visiting a particular basin under uniform sampling is proportional to the energy of that minima. Therefore, the hypervolume of each energy basin is proportional to the basin depth. This fact can be derived by noting that Eqn. (2) is a quadratic form, so the Hessian matrix is a constant. Therefore, all local minima have the same mean curvature. Given that low lying energy minima are wide, we can extract the position of energy minima in the space of input variables by studying the correlation structure of the binary dataset. We note that the correlation between basin hypervolume and basin depth appears in many complex physical systems beyond the Ising model Doye et al. (1998); Doye and Massen (2005); Pickard and Needs (2011).

A lingering question is whether our inference procedure is robust to the choice of threshold. To test this, we consider $m$ Hopfield patterns, chosen as eigenvectors of a symmetrized $p\times p$ Gaussian random matrix, with the Hopfield energy chosen to be Gaussian distributed with mean $10$ and unit variance. We draw 10000 samples randomly and compute the correlation matrix with the lowest energy $N$ samples. Figure 3 shows that our method is robust: the correlation coefficient between $J$ and $J_{\mathrm{MP}}$ is large and constant for a wide range of thresholds and number of Hopfield patterns. The question of how many energy minima can be recovered from the binary data and a thermodynamic interpretation is discussed in the Supplemental Information.

We now turn to consider two common scenarios that break the assumptions made so far, stratified sampling and more complex energy landscapes.

Stratified sampling: Thus far we assumed that the sampling is uniform before thresholding. However, in many settings, the sampling is biased. To model this effect, we draw 5000 random samples, but freeze the first 5 variables to $+1$ for the first 2500 samples and the last 5 variables to $-1$ for the remaining 2500 samples. We then evaluate the energy, and take again the lowest $10^{\mathrm{th}}$ percentile. Figure 4A-D shows that the frozen variables introduce sample-sample correlations, and now there are 4 significant eigenvectors with the first Hopfield pattern demoted to the second largest eigenvector (Figure 4C). As such, the informative eigenvectors are still present in $J_{\mathrm{MP}}$, but the eigenvalues are misplaced.

In this case, naïve random matrix cleaning does not recover $J$ (Figure 4E), since there is no a priori reason to discard the first eigenvector unless we know the Hopfield patterns beforehand. We need additional information – for which we turn to the quantitative measurements – to accurately recover the Hopfield energies. Figure 4F shows that an additional 500 quantitative measurement allow us to recover the coupling matrix (MAE = 0.12) using ridge regression and the ansatz Eqn. (3). The error is significantly larger (MAE = 0.31) if only the quantitative measurements are used.

Complex energy landscapes: The geometric property that the depth of an energy minima is related to its hypervolume is not universal to all energy landscapes Wales et al. (1998); Wales (2003). A natural question is whether the significant eigenvectors and eigenvalues of the correlation matrices of samples below/above an energy threshold allow us to infer features of a complex energy landscapes. We consider a landscape that comprises a sum of Gaussians

[TABLE]

This landscape has the property that the depth of each energy minima, $E_{i}$ (located at $\boldsymbol{\zeta}_{i}$), is inversely proportional to its width $1/E_{i}$. As above, we let $(E_{1},E_{2},E_{3})=(-30,-25,-20)$ and generate Hopfield patterns by diagonalising a symmetrized Gaussian random matrix. We draw $5000$ samples and threshold to find the $500$ lowest energy samples. Figure 5 shows that there are again three significant eigenvectors above the Marčenko-Pastur threshold, but the lowest energy Hopfield pattern is demoted to the third eigenvector, while the highest energy Hopfield pattern is promoted to the top eigenvector. This is expected: the eigenvalue corresponding to each minimum is proportional to the number of samples near that minimum, i.e. the basin volume, which in this case is not proportional to basin depth. However, the eigenvectors still indicate the locations of the energy minima, motivating the approach described in Eqn. (2), where we use these eigenvectors to identify features.

In conclusion, we develop a general strategy, grounded in statistical physics, which integrates coarse and fine measurements to yield a predictive model. Since coarse measurements are often significantly less costly to obtain, our strategy provides a new avenue for experiment design. Although our Letter only considered an Ising-type model, the fact that the eigenvectors of the correlation matrix of coarse measurements point toward energy minima suggests a natural way to integrate our result into more complex non-linear models, for example by using $\mathbf{f}\cdot\mathbf{u}_{i}$, the overlaps between the sample vector and each eigenvector, as inputs to a general nonlinear function such as an artificial neural network.