The problem of spurious correlations is fundamental in statistics. Spurious correlations may be due to confounding or selection bias. Confounding is induced by a common unobserved cause that underlies both the predictor variables and the traits. Selection bias emerges arises from taking non-random subsets of the population, where some members are less likely to be subsampled than others [ 19 ]. Both e ects may result in the phenomena of spurious correlations which we can treat here simultaneously [ 5, 6, 7 ]. Population structure [ 9 ] implies that due to common ancestry, genes of individuals that are related co-inherit a large number of genes, making them more similar to each other, whereas the genes of people of unrelated ancestry are obtained independently, making them more dissimilar. Population structure is the root of many unwanted biases. For example, when data is collected only in selected geographical areas (such as in specific hospitals), one thereby introduces a selection bias into the sample, meaning that the collected genes do not represent the overall population. This can heavily distort the prediction quality of a classifier [ 13 ]. Another problem is that people who live geographically close often share other factors, such as similar environmental factors or culture. This, in turn, can lead to similar phenotypes (such as overweight, drinking habits, or diabetes). Thus, because genes correlate with location and location may correlate with specific phenotypes, there is a resulting correlation between genes and these phenotypes that does not have a causal interpretation—another manifestation of confounding by population structure [ 20 ]. It is an active area of research to find models that are less prone to spurious correlations [ 13 ]. In this paper, we present such a model for the setup of binary classification.

A popular approach to correcting for spurious correlations relies on similarity kernels, or kinship matrices [ 9 ]. Given n samples, we can construct an n n matrix K that quantifies the similarity between samples based on some arbitrary measure. In the case of confounding by population structure, one typically chooses Ki j = Xi>X j, where Xi is a vector of genetic features of individual i. As Ki j contains the scalar products between the genetic vectors of individuals i and j, it is a sensible measure of genetic similarity. As another example, when correcting for confounding by age, then we can choose K to be a matrix that contains 1 if two individuals have the same age, and zero otherwise. Details of constructing similarity kernels can be found in [ 9 ]. Next, we explain how the similarity matrix can be used to correct for confounding. 2.2

Our model builds on the LMM-Lasso [ 18 ], an important method of statistical genetics to limit the impact of confounding. While the LMM-Lasso relies on linear regression, we generalize this approach to the much more involved classification setup, where the target values are binary. The correlated probit model is yi = sign Xi>w + i ; = ( 1; : : : ; n)>

N (0; ): (1) In the special case of = I, this is just the ordinary (uncorrelated) probit model. In the following, we refer to this model as Correlated Probit Regression (CPR). For now, we assume that the covariance matrix is fixed and known. In our empirical studies we use the parametrization Eq. 13 where the parameters are estimated from the data. We now derive an objective function to find an estimate of the model parameter w. To simplify the notation, we will without loss of generality assume that all observed binary labels yi are 1. The reason why this assumption is no constraint is that we can always perform a linear transformation to absorb the sign of the labels into the data matrix and noise covariance1. Thus, when working with this transformed data matrix and noise covariance, our assumption is satisfied.

The likelihood function, thus, is the probability that all transformed labels are 1. This is satisfied when Xi>w + i > 0. When integrating over all realizations of noise, the resulting (marginal) likelihood is

P(8i : yi = 1jw) = =

P(8i : Xi>w + i > 0jw) Z

Rn N ( ; X>w; ) dn : + (2) The marginal likelihood is hence an integral of the multivariate Gaussian over the positive orthant. In Section 3, we will present e cient approximations of this integral. Before we get there, we further characterize the model. Next, we turn the correlated probit model into a model for feature selection. We are interested in a point estimate of the weight vector w that is sparse, i.e. contains zeros almost everywhere. This is well motivated in statistical genetics for phenotypes or diseases that are believed to be caused by a small number of genes. Sparsity is achieved using the Lasso [ 15 ], where we add an `1-norm regularizer to the negative marginal likelihood: The fact that the noise variable and the weight vector w have di erent priors or regularizations makes the model identifiable and lets us cleanly distinguish between linear e ects and e ects of correlated noise. It is easy to show that the objective function Eq. 3 is convex. 2.3

Before we come to inference, we point out how our approach connects to other methods. When removing the probit likelihood, the model becomes the LMM-Lasso [ 18 ], hence P(Yjw) = N (Y; X>w; ). This model has shown to improve selection of true non-zero e ects as well as prediction quality [ 18 ]. Our model is a natural extension of the LMM-Lasso to binary outcomes, such as the disease status of a patient. As we explain in this paper, inference of our model is, however, much more challenging than in [ 18 ]. 1To this end, we apply the transformations X diag(y) diag(y). diag(y)X and

Furthermore, by construction, our model captures two limiting cases: uncorrelated probit regression and Gaussian process (GP) classification. To obtain uncorrelated probit regression, we simply assume a covariance matrix proportional to unity. To obtain GP classification, we simply omit the fixed e ect (i.e., we set w = 0 in Eq. 2) so that our model likelihood becomes P(Y = Yobsjw) = RRn+ N ( ; 0; ) dn , where hence the noise variable plays the role of the latent function f in GPs [ 16 ]. We will compare our method to all three related models in the experimental part of the paper and show enhanced accuracy.

CPR is a EM-type algorithm [ 27 ]. In the outer loop (the Mstep), we follow gradients to optimize the objective. Since this objective function has an `1-norm regularizer, we have to split this outer optimization routine into two parts, one that optimizes the likelihood and one that optimizes the regularizer. This is described in 3.1.1. The inner loop (the E-step) consists of computing the gradient and the Hessian of the likelihood term by means of approximate inference, which is described in Sections 3.1.2 and 3.1.3. 3.1.1 The `1-norm in the objective function Eq. 3 prevents us from directly applying gradient based methods such as Newton’s method. A solution is given by ADMM that involves a generalized objective:

L(w; z; ) := log

Z

RN N ( ; X>w; )dn + 0jjzjj11 + + >(w

1 z) + 2 cjjw zjj22: We minimize over w and z and maximize over . In alternating between the minimization updates for w, z and a gradient step in , we solve the original problem [ 22 ]. While the updates for z and have analytic solutions, we compute the updates for w by numerical optimization. The part of the ADMM objective L(w; z; ) depending on w, called L(w) for brevity, is e ectively `2-norm regularized, enabling us to compute the gradient and the Hessian. This allows us to apply Newton’s Method to obtain the ADMM update in w. 3.1.2

The inner loop of the EM-algorithm amounts to computing the gradient and Hessian of L(w; z; ). These are not available in closed-form, but in terms of the first and second moment of a truncated Gaussian density.

Since computing the derivates of the linear and quadratic term is straightforward, we focus on L0(w) := log RR+N N ( ; X>w; )dn , which contains the intractable integral. In the following, we use the short hand notation (w) = X>w: (4) (5) (6) means of the posterior (the truncated Gaussian) and the untruncated Gaussian. The gradient and Hessian of L0(w) are given by rwL0(w) =

H0(w) =

X[ 1X>; 1( p >) 1 (7) Note that the variable w enters through p(w) and (w). Next, we describe how we approximate the intractable expectations involved in Eq. 7. 3.1.3

In Eq. 7 we have expressed the gradient and Hessian of L0(w) in terms of the first and second moment of the posterior Eq. 5. The problem is that computing the moments involves intractable expectations over this distribution. We employ EP [ 17 ] to approximate these expectations. Note that also other approximate inference schemes are possible, such as variational inference or sampling methods [ 16 ]. EP approximates moments of the posterior p( j ; ) in terms of a variational distribution q( ), approximately minimizing the Kullback-Leibler divergence, q ( j q ; q ) = arg min Ep[log p( j ; )] q

Ep[log q( j q; q)] : The variational distribution q ( ) is an un-truncated Gaussian q ( ; q ; q ) = N ( ; q ; q ), characterized by the variational parameters q and q . We approximate the mean and covariance of the posterior p in terms of the variational distribution, p q , and p q . We warm-start each gradient computation with the optimal parameters of the earlier iteration.

Algorithm 1 summarizes our procedure. We denote the EP algorithm for approximating the first and second moment of the truncated Gaussian by EP( ; ). Here, and are the mean and covariance matrix of the un-truncated Gaussian. The subroutine returns the first and second moments of the truncated distributions q and q. When initialized with the outcomes of earlier iterations, this subroutine converges within a single EP loop.

Our algorithm thus consists of three nested loops; the outer ADMM loop, containing the Newton optimization loop for computing the update in w and the inner EP loop, which computes the moments of the posterior. We choose stopping criterion 1 to be the convergence criterion proposed by Boyd [ 22 ] and choose criterion 2 to be always fulfilled, i. e. we perform only one Newton optimization step in the inner loop. Our experiments showed that doing only one Newton optimization step, instead of executing until convergence, is stable and leads to great speed improvements. ADMM is known to converge even when the minimizations in the ADMM scheme are not carried out exactly (see e.g. [ 28 ]).

It is convenient to introduce the following probability distribution: p( j ; ) = 1[ 2 Rn+] N ( ; ; )

RRn+ N ( ; ; ) dn Above, 1[ ] is the indicator function. Eq. 5 is just the multivariate Gaussian, truncated and normalized to the positive orthant; we call it the posterior distribution. We furthermore introduce p(w) = p(w) =

Ep( j (w); ) [ ] ; Ep( j (w); ) h( p(w))( p(w))>i : p and p are the mean and the covariance of the truncated multivariate Gaussian, as opposed to and which are the non-truncated ones.

In the following we abbreviate p p(w) and p p(w), and write = p for the di erence between the

X = y X˜ nnpre-process the data

= diag(y) ˜ diag(y) repeat nnget wk+1 by EP and Newton’s Method initialize w = wk repeat ( q; q) EP(X>w; )

= q X>w g = > 1 X> + c(w zk + k)> H = X[ 1 1( q >) w = w t H 1g until criterion 2 nnADMM updates wk+1 = w zk+1 = S =c(wk+1 + k) nnsoft thresholding k+1 = k + wk+1 zk+1 until criterion 1 1]X> + cI 3.2

For simplicity, we use covariance matrices of a special structure, which allows us to derive an alternative formulation of the correlated probit model2. In particular, we assume that is a combination of diagonal noise and a linear kernel of the data matrix, =

1I + 2 X> X: The linear kernel X> X measures similarities between genes and therefore models the e ect of genetic similarity between samples due to population structure. We use the following Gaussian integral identity: (8) (9) marginalized out), it correlates the samples. We interpret w0 as a confounder.

The simplest approximation to Eq. 9 is to substitute the integral over w0 by its maximum a posteriori (MAP) value, leading to the new objective function:

L(w; w0) = + n X log Under the MAP approximation, the likelihood contribution to the objective function becomes completely symmetric in w and w0: only the sum w + w0 enters. The di erence between the two weight vectors w and w0 in this approximation is only due to the di erent regularizers: while w0 has an `2-nrom regularizer and is therefore dense, w is `1-norm regularized and therefore sparse.

For optimizing the MAP approximated objective function Eq. 10 jointly in w and w0, we introduce a block coordinate descent scheme alternating between updates in w and w0. For updating w0 we use gradient descent, while for updating w we employ ADMM (c.f. Section 3.1.1). Note that the procedure could be made faster by using a second-order optimization method for obtaining the updates in w0. Under the MAP approximation, every feature gets a small non-zero weight from w0, and only selected features get a stronger weight from w. The idea is that w0 models the population structure, which a ects all genes. In contrast, we are interested in learning the sparse weight vector w, which has a causal interpretation because it involves only a small number of features.3 The MAP approximation is computationally more convenient, but it has its limits. In the original correlated probit model in Eq. 1, we marginalize over the confounder, which is more expensive. In contrast, under the MAP approximation we optimize over w0 and the the objective function factorizes over n, which means that we have broken the correlations between the samples. This comes at the cost of reduced prediction performance. Since the MAP estimate of the confounder does not capture the full information of its distribution, the MAP probit model tends to generalize not as well as the (full) correlated probit model. We compare both approaches experimentally in Section 4. 3.3

Stochastic gradient Monte Carlo methods are an active area of research in scalable Bayesian inference. These methods approximately sample from a posterior by using only a 3Note that the interplay of two weight vectors is di erent from an elastic net regularizer. Above, ( ) is the cumulative standard normal distribution function. We have introduced the new Gaussian noise variable w0. Conditioned on w0, the remaining integrals factorize over n. However, since w0 is unobserved (hence 2Note that the approach can be easily generalized to arbitrary covariance matrix by considering the Cholesky decomposition = BB>. subset of data for generating a sample and, therefore, being scalable to big datasets. This is done by using stochastic optimization to provide e cient proposals for MetropolisHastings algorithms with high acceptance rates. We propose two versions of SG-CPR, one builds on Stochastic Gradient Langevin Dynamics (SGLD) [ 24 ] and the other on Constant Stochastic Gradient Descent (c-SGD) [ 25 ]. These methods assume that the likelihood factorizes, conditioned on the global variables. Up to a constant, the log posterior is

n log p( jy) =c X log p(yij ) + log p( ):

i=1 Let S be a set of S random indices drawn uniformly at random from the index set f1; : : : ; ng. The stochastic gradient w.r.t. the minibatch S of the log likelihood term is gˆS( ) =

S1 X r log p(yij ):

i2S SGLD and c-SGD work as follows. SGLD performs decreasing stochastic gradient step on the negative log joint distribution, but adds artificial noise to prevent convergence to the optimum. Instead, the algorithm converges to a stationary distribution, which can be shown to be the posterior [ 24 ]. Constant SGD, on the other hand, only approximates the posterior. It converges faster because it operates with constant step sizes.

The following formula summarizes the two methods: SGLD : c-SGD: i+1 = i+1 = i + i + 2t gˆS( t) + t; gˆS( t); t

N (0; t); (11) where t is a suitable decreasing learning rate for SGLD. For c-SGD, it was shown in [ 25 ] that the optimal constant learning rate that best approximates the posterior equals 2dS nTr(E[gˆgˆT ]) . This definition involves the minibatch size S , feature space dimension d, and the stochastic gradient noise covariance near the optimum, E hgˆgˆT i.

We now explain how SGLD and c-SGD can be used for inference in the correlated probit model. Recall that the aim is to find the MAP estimate of the model by optimizing the objective function L(w) Eq. 3. In Section 3.2, we introduced an auxiliary variable w0 and obtained the identity c L(w) = log

Z

Rd which involves the posterior of the confounder p(w0jy; w). This posterior can be approximated using SGLD or c-SGD. Since this algorithm is based on stochastic gradient descent SG-CPR scales to hundreds of thousands data points.

For the real-world and synthetic experiments, we first need to make a choice for the class of kernels that we use for the covariance matrix. We choose a combination of three contributions, = 1I + 2 X>X + 3 side: (13) The third term is optional and depends on the context; it is a kernel that we extract from side information in the form of an additional feature matrix X0, where we choose side as an RBF kernel [ 29 ] on top of the side information X0. Note that this way, the data matrix enters the model both through the linear e ect but also through the linear kernel. We evaluate the methods by using n instances of the dataset for training and splitting the remaining examples equally into validation and test sets. This process is repeated 50 times, over which we report on average accuracies or areas under the Receiver Operating Characteristic (ROC) curve (AUCs) as well as standard errors [ 30 ].

The hyper parameters k in Eq. 13, together with the regularization parameter 0, were determined on the validation set, using grid search over a su ciently large parameter space (optimal values are attained inside the grid; in most cases k 2 [0:1; 1000] for k = 0; 1; 2; 3). For all datasets, the features were centered and scaled to unit standard deviation, except in Section 4.4, where they are binary. In Sections 4.3 and 4.4, we show that including a linear kernel into the covariance matrix leads to top features which are less correlated with the population structure in comparison to the features of uncorrelated probit regression. The correlation plots4 in Fig. 4 show the mean correlation of the top features with population structure and the corresponding standard errors. 4.2

We generated n = 200 synthetic data points in d = 50 dimensions as follows. We generate a weight vector w 2 Rd with k entries being 1, and the other d k entries being 0, where 1 k d. We then create a random covariance matrix side 2 Rn n, which serves as side information matrix5. We draw n points X = fx1; : : : xng independently from a uniform distribution over the unit cube [ 1; 1]d and create the labels according to the probit model Eq. 1, using side as covariance matrix. We reserve 100 samples for training and 50 for validation and testing, respectively. As a benchmark we introduce the oracle classifier, where we use the correlated probit model (with covariance matrix side) but skip the training and instead use the true underlying w for prediction.

For several 1 k d, we generate a dataset according to the above described procedure. In Fig. 1, we report on the so-achieved accuracies with respect to the percentage of non-zero features ( dk ). We observe that in the sparse scenarios ( 20% non-zero features), GP classification and LMM-Lasso are clearly outperformed by CPR, achieving an accuracy up to 10 percentage points and 23 percentage points higher, respectively. Due to being `1-norm regularized and therefore, having the capability of exploiting sparsity, uncorrelated probit regression performs best in this regime among the competitors, but still substantially worse than CPR. LMM-Lasso is also `1-norm regularized but is not designed for a classification setting. Therefore, it cannot beat uncorrelated probit regression. In the dense scenarios, CPR outperforms LMM-Lasso (by 1 to 4 percentage points) and performs similarly well to GP classification, which also takes the correlation structure into account. In this scenario, uncorrelated probit regression is 4The correlation plots in Fig. 4 are created according to [ 8 ] as follows. First, we randomly choose 70% of the available data as training set and obtain a weight vector w by training. We compute the empirical Pearson correlation coe cient of each feature with the first principle component of the linear kernel on top of the data. This is a way to measure the correlation with the population structure [ 31 ]. We define the index set I by taking the absolute value of each entry of w and sorting them in descending order. We now sort the so-obtained list of correlation coe cients with respect to the index set I and obtain a resorted list of correlation coe cients (c1; : : : ; cn). In the last step, we obtain a new list (cˆ1; : : : ; cˆn) by smoothing the values, computing cˆi := 1i Pik ck. Finally, we plot the values (cˆ1; : : : ; cˆn) with respect to I. This procedure was repeated 30 times for di erent random choices of training sets.

5The covariance matrix was created as follows. The random generator in MATLAB version 8.3.0.532 was initialized to seed = 20 using the rng(20) command. The matrix side was realized in two steps via A=2*rand(50,n)-1 and side=3*A'*A+0.6*eye(200)+3*ones(200,200). clearly worse than the other methods, because it does not take the correlation structure into account. We observe that in all scenarios the prediction performance of CPR-MAP is between uncorrelated probit regression and CPR. In Fig. 2, we inspect the computed feature weights (green dots) of `1-norm regularized and `2-norm regularized CPR, respectively. The blue solid line represents the ground truth (the true underlying weight vector w with k = 10 entries non-zero). We observe that the `1-norm regularized probit model finds the true weights without su ering from large noise as the `2-norm regularized counterpart does.

Toy Experiment

Oracle CPR CPR-MAP Uncorrelated Probit GP

LMM-Lasso 20%

40% 60% percentage of non-zero features 80% 100% We obtained the dataset by [ 32 ] from the National Center for Biotechnology Information website 6, which includes 40 blood samples from patients with active tuberculosis as well as 103 healthy controls, together with the transcriptional signature of blood samples measured in a microarray experiment with 48,803 gene expression levels, which serve as features for our purposes. Also available is the age of the subjects when the blood sample was taken, from which we compute side7. All competing methods are trained by using various training set sizes n 2 [40; 80]. To be consistent with previous studies (e. g. [ 8 ]), we report on the area under the ROC curve (AUC), rather than accuracy, where we vary the hyperparameters k. The results are shown in Fig. 3, left.

We observe that CPR achieves a consistent improvement over its uncorrelated counterpart (by up to 12 percentage 6http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE19491

7We compute side as RBF kernel on top of the side information age using bandwidth = 0:2. points), GP classification (by up to 3 percentage points), LMM-Lasso (by up to 7 percentage points) and CPR-MAP (by up to 7 percentage points). In Fig. 4, left, we show the correlation of the top features with population structure (as confounding factor) for correlated and uncorrelated probit regression. The plot was created as explained in Section 4.1. We find that the features obtained by CPR show much less correlation with population structure than the features of uncorrelated probit regression. By inspecting the correlation coe cients of the first top 10 features of both methods, we observe that the features found by CPR are 40 % less correlated with the confounder. This is because population structure was built into our model as a source of correlated noise.

AUCs in the Tuberculosis Dataset CPR CPR-MAP Uncorrelated Probit GP

LMM-LASSO 6540 45 50

55 60 65 number of training samples 70 75 80

ROC Curves for the Malware Dataset 90 85 80 C U A 75 70 1.0 0.8 e tra0.6 e v iit s o p -e0.4 u tr 0.2 CPR CPR-MAP Uncorrelated Probit GP LMM-Lasso 0.4 0.5 0.00.0 0.1 0.2 0.3 false-positive rate We experiment on the Drebin dataset8 [ 33 ], which contains 5,560 Android software applications from 179 di erent malware families. There are 545,333 binary features; each feature denotes the presence or absence of a certain source code string (such as a permission, an API call or a network address). It makes sense to look for sparse feature vector [ 33 ], as only a small number of strings are truly characteristic of a malware. The idea is that we consider 8http://user.informatik.uni-goettingen.de/ ~darp/drebin/download.html 20 40 60 80 index sorted by decreasing absolute feature weight

Malware Dataset

CPR

Uncorrelated probit 20 40 60 80 index sorted by decreasing absolute feature weight populations of di erent families of malware when training, and hence correct for the analogue of genetic population structure in this new context, that we call “malware structure”. We concentrate on the top 10 most frequently occurring malware families in the dataset9. We took 10 instances from each family, forming together a malicious set of 100 and a benign set of another 100 instances (i.e., in total 200 samples). We employ n = 80 instances for training and stratify in the sense that we make sure that each training/validation/test set contains 50% benign samples and an equal amount of malware instances from each family. Since no side information is available, we only use a linear kernel and the identity matrix as components for the correlation matrix. We report on the (normalized) area under the ROC curve over the interval [0; 0:1] and denote this performance measure by AUC0:1. In Fig. 3, right, we show the ROC curves and in Table 1 the achieved AUC0:1. We observe that correlated probit regression (CPR) achieves a consistent improvement in terms of AUC0:1 over the competitors (by up to 8.4 percentage points). Furthermore, in Fig. 4, right, we plot the correlation of the top features of correlated and uncorrelated probit regression with population structure. We observe that CPR leads to features, which are much less correlated with the malware 9Geinimi, FakeDoc, Kmin, Iconosys, BaseBridge, GinMaster, Opfake, Plankton, FakeInstaller, DroidKungFu. 4.5

We experiment on genotype and phenotype data consisting of 199 genetically di erent samples from the model plant Arabidopsis thaliana [ 34 ]. The genotype of each sample comprises 216,130 single nucleotide polymorphism (SNP) features. The phenotype that we aim to predict is early or late flowering of a plant when grown at ten degrees centigrade. The original dataset contains the flowering time for each of the 199 genotypes. We split the dataset into the lower and upper 45%-quantiles of the flowering time and removed the middle 10%. We then binarized the labels, resulting in a set of 180 instances from which we use n = 150 instances for training. The results are reported in Table 2

CPR 84:1 0:2

CPR-MAP 83:6 0:3

Uncorr. Probit 83:5 0:2

GP 83:6 0:2

LMM-Lasso 79:7 0:2 and show that CPR has a slight advantage of at least 0:5 percentage points in AUC over the competitors. An analysis restricted to the ten SNPs with largest absolute regression weights in our model showed that they lie within four well-annotated genes that all convincingly can be related to flowering, structure and growth: the gene AT2G21930 is a growth protein that is expressed during flowering, AT4G27360 is involved in microtubule motor activity, AT3G48320 is a membrane protein, involved in plant structure, and AT5G28040 is a DNA binding protein that is expressed during flowering. 5

We thank Manfred Opper, Mehryar Mohri, David Blei, Rajesh Ranganath, Maja Rudolph, and Gunnar Ra¨tsch for stimulating discussions. SM acknowledges the support of the U.S. National Science Foundation I2CAM International Materials Institute Award, Grant DMR-0844115, and the NSF Schloss Dagstuhl support grant for junior researchers (CNS-1257011). SM and MK gratefully acknowledge the support of the NVIDIA Corporation for the donation of the Tesla K40 GPU. MK acknowledges support from the German Research Foundation (DFG) award KL 2698/2-1 and from the Federal Ministry of Science and Education (BMBF) award 031L0023A. preprint arXiv:1106.1157, 2011. [39] M. W. Seeger and H. Nickisch, “Large scale bayesian inference and experimental design for sparse linear models,” SIAM Journal on Imaging Sciences, vol. 4, no. 1, pp. 166– 199, 2011.