{'fraction_non_alphanumeric': 0.04653715289210844, 'fraction_numerical': 0.014103830534518714, 'mean_word_length': 4.282649659370923, 'pattern_counts': {'":': 8, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-Sample Testing in High-Dimensional Models

We propose novel methodology for testing equality of model parameters between two high-dimensional populations. The technique is very general and applicable to a wide range of models. The method is based on sample splitting: the data is split into two parts; on the first part we reduce the dimensionality of the model to a manageable size; on the second part we perform significance testing (p-value calculation) based on a restricted likelihood ratio statistic. Assuming that both populations arise from the same distribution, we show that the restricted likelihood ratio statistic is asymptotically distributed as a weighted sum of chi-squares with weights which can be efficiently estimated from the data. In high-dimensional problems, a single data split can result in a"p-value lottery". To ameliorate this effect, we iterate the splitting process and aggregate the resulting p-values. This multi-split approach provides improved p-values. We illustrate the use of our general approach in two-sample comparisons of high-dimensional regression models ("differential regression") and graphical models ("differential network"). In both cases we show results on simulated data as well as real data from recent, high-throughput cancer studies.


Introduction and Motivation
We consider the general two-sample testing problem where the goal is to test whether or not two independent populations U and V , parameterized by φ u and φ v respectively, with φ u , φ v ∈ Φ ⊂ R p , arise from the same distribution. The hypothesis testing problem of interest is (1.1) In this paper the focus is on the high-dimensional setting where the number of available samples per population (n u and n v ) is small compared to the dimensionality of the parameter space p.
In a setup where n u and n v are much larger than p the ordinary likelihood-ratio test offers a very general solution to problem (1.1). It is well-known that under H 0 the likelihood ratio statistic is asymptotically χ 2 p distributed which then allows computation of confidence intervals and p-values. However, in settings where p is large compared to sample sizes, the likelihood-ratio statistic is ill-behaved and the classical asymptotic set-up cannot be relied upon to test statistical significance.
Our approach for solving the general high-dimensional two-sample problem (1.1) is motivated by the screen and clean procedure (Wasserman and Roeder, 2009) originally developed for improved variable selection in the high-dimensional regression model. The idea is to split the data from both populations into two parts; perform dimensionality reduction on one part and pursue significance testing on the other part of the data. In more detail, on the first split we select three subsets of Φ by screening for the most relevant parameters of population U and V individually but also by selecting the important parameters of the pooled data from both populations. In this way we obtain an individual model with different parameter spaces for U and V and a joint model which describes both populations by a common parameter space. On the second split of the data, we then can compute the likelihood ratio statistic between these two models, the so-called restricted likelihood ratio statistic. A crucial observation is that the two models are non-nested and therefore nonstandard tools are needed to obtain the asymptotic null distribution. We apply the theory on model selection and non-nested hypotheses developed by Vuong (1989) to the twosample scenario and show that the null distribution asymptotically approaches a weighted sum of independent chi-squared random variables with weights which can be efficiently estimated from the second split of the data. Importantly, the weighted sum of chi-squared approximation is invoked only in the second split, where the model dimensionalty has been reduced.
As indicated above our approach involves a screening or model selection step prior to significance testing. Our method requires that the parameter set selected by a specific screening procedure contains the true model parameter (screening property) and that the selected set is small compared to the sample size (sparsity property). The first property is needed for deriving the asymptotic null distribution, whereas the latter property justifies asymptotic approximation. The two conditions which we impose here are much weaker than requiring consistency in variable selection. We propose to use ℓ 1 -penalized maximum likelihood estimation (Tibshirani, 1996;Fan and Li, 2001) with regularization parameter chosen by cross-validation. This sets automatically non-relevant parameter components to zero and allows for efficient screening even in scenarios with very large p. ℓ 1 penalization leads to sparse estimates and there is theoretical evidence that the screening property holds under mild additional assumptions.
In current applied statistics, a wide range of applications are faced with high-dimensional data. Parameter estimation in the "large p, small n" setting has been extensively studied in theory (Bühlmann and van de Geer, 2011) and also applied with success in many areas. Since interpretation of parameters is crucial in applied science there is a clear need to assess uncertainty and statistical significance in such settings. However, significance testing in high-dimensional settings has only recently attracted attention. Meinshausen et al. (2009) and Bühlmann (2012) consider testing in the high-dimensional linear model. In the context of high-dimensional two-sample comparison Bai and Saranadasa (1996), Chen and Qin (2010), Lopes et al. (2012) and others treat testing for differences in the population means. And recently, Cai et al. (2011) and Li and Chen (2012) developed a two-sample test for covariance matrices of two high-dimensional populations.
The approach proposed in this paper tackles the high-dimensional two-sample testing in a very general setting. In our methodology φ u and φ v can parameterize any model of interest. In the empirical examples we show below we focus on two specific applications of our approach: (i) High-dimensional regression, where two populations may differ with respect to regression models and (ii) Graphical modelling, where the two populations may differ with respect to conditional independence structure. In analogy to the term "differential expression" as widely-used for testing means in gene expression studies, we call these "differential regression" and "differential network" respectively. Both highdimensional regression and graphical models are now widely used in biological applications, and very often scientific interest focuses on potential differences between populations (such as disease types, cell types, environmental conditions etc.). However, to date in the highdimensional setting, two sample testing concerning such models has not been well studied. The methodology we propose offers a way to directly test hypotheses concerning differences in molecular influences or biological network structure using high-throughput data.
Note that x can be empty (i.e., l = 0). Let population U = (Y u , X u ) and V = (Y v , X v ), where X u , X v are both generated from the same distribution and conditional on the X's, Y u and Y v are generated from d(·|x; φ u ) and d(·|x; φ v ) respectively.
The goal is to solve the general two-sample testing problem (1.1) given data-matrices U = (Y u , X u ) ∈ R nu×(k+l) and V = (Y v , X v ) ∈ R nv×(k+l) , representing i.i.d. random samples of U and V . We consider the high-dimensional case with p ≫ min{n v , n u } and assume that the data generating parameters φ u and φ v are sparse, meaning that many of their components are equal to zero. If k = 1, l = 0, we have a classical (univariate) two-sample set-up. If k = 1 and d(y|x; φ) is the univariate Normal distribution with mean β T x and noise variance σ 2 , φ = (β, σ 2 ), then U and V follow linear regression models. If l = 0 and d(y; φ) is the multivariate Normal distribution with φ representing the inverse covariance matrix, then U and V follow a Gaussian graphical model and (1.1) asks whether or not two graphical models (in short: networks) are significantly different. We will treat these two examples in detail in Section 3 and 4. We refer to the regression case as differential regression and to the graphical model case as differential network.
Our methodology for solving (1.1) is based on sample splitting and has its inspiration from the work by Wasserman and Roeder (2009). We randomly divide data from population U into two parts U in and U out of equal size and proceed in the same way with population V which yields V in and V out . In a first step the dimensionality of the parameter space Φ is reduced by filtering out (potentially many) redundant components. We do this by applying a screening procedure I on U in and V in separately, but also on pooled data (U in , V in ). In this way we obtain models with lower dimensionality than the full dimension p. The first model describes populations U and V individually using different reduced-parameter spaces, whereas the second model explains both populations jointly with a single reducedparameter space. In a second step, we then evaluate the restricted likelihood-ratio statistic on the held-out data U out and V out and perform significance testing.

Screening and ℓ 1 -regularization
We have p ≫ n and φ o is supposed to be sparse.
A screening procedure I selects, based on data Z, a set A of active parameter components by a map I(Z) ⊂ {1, . . . , p}. The map I then defines the active parameter space by There are two basic requirements on the screening procedure I. Firstly, the procedure should get rid of many non-relevant components in the sense that the set of active parameter components I(Z) should be small compared to n. Secondly, the active parameter space Φ I(Z) should contain the true parameter φ o . We refer to these requirements as the sparsity and screening property: • Sparsity property: |I(Z)| is small compared to n.
• Screening property: Formulating the sparsity property precisely, i.e. specifying the rate at which the size of the active set I(Z) can grow with n, is a research topic in its own right (see also the discussion in Section 5). We do not address this topic here. The screening property guarantees that the model {d(·|x; φ) : φ ∈ Φ I(Z) }, selected by I, is correctly specified in the sense that it contains the true density function d(·|x; φ o ). L1-penalized likelihood methods (Fan and Li, 2001) serve as our prime example of screening procedures. Consider estimators of the form where ℓ(φ, Y|X) denotes the conditional log-likelihood and λ is a non-negative regularization parameter. In the context of linear regression (2.3) coincides with the Lasso estimator introduced by Tibshirani (1996). Other important example of the form (2.3) include the elastic net (Zou and Hastie, 2005), the grouped Lasso (Yuan and Lin, 2006), the graphical Lasso (Friedman et al., 2008) or the Lasso for generalized linear models (Park and Hastie, 2007). Estimators of the form (2.3) have been shown to be extremely powerful: they are suitable for high-dimensional data and lead to sparse solutions. A screening procedure can be defined by setting The question is whether or not I λ (Z) satisfies both the sparsity and the screening property.
In principle the answer depends on the choice of the tuning parameter λ. Too strong regularization (too large λ) leads to very sparse solutions, but with relevant parameter components incorrectly set to zero, whereas too little regularization (too small λ) results in too large active sets. Common practice is to choose λ in a prediction optimal manner, for example by optimizing a cross-validation score. This gives typically sparse solutions with a larger number of non-zero components than the true number (Meinshausen and Bühlmann, 2006). It seems to be reasonable in practice to assume that the sparsity and screening property hold for the procedure I λ (·) obtained via ℓ 1regularization. In fact, Bühlmann (2012) points out that the screening property is a very useful concept for ℓ 1 -regularization which holds under much milder assumptions than consistency in variable selection.
By applying the screening procedure I λ on U in and V in separately, we obtain active-sets Further we obtain active-set by applying I jointly on (U in , V in ). For the rest of this section we treat the active-sets I u , I v and I uv as fixed, i.e., not depending on the sample size.

Restricted likelihood ratio statistic
We define model M joint with shared parameter space Φ Iuv for both population U and V jointly as and model M ind with individual parameter spaces Φ Iu and Φ Iv for each population as The log-likelihood functions with respect to models M joint and M ind are given by We will test hypothesis (1.1) based on the restricted log-likelihood ratio test-statistic defined as where (φ u ,φ v ) andφ uv denote the maximum likelihood estimators corresponding to models M ind and M joint .
It is crucial to note that testing based on (2.5) is non-trivial and involves non-nested model comparison (Vuong, 1989). The difficulty arises from the fact that model M joint is not nested in M ind , or in other terms Φ J ⊂ Φ Iu ∪ Φ Iv . Non-nestedness of these models is fundamental and is not only an artefact of the random nature of the screening procedure.
If φ u = φ v , than I applied on two different populations mixed together will set different components of φ to zero than I applied two both populations individually. This is a consequence of model misspecification and is also well known as Simpson's paradox in which an association present in two different groups can be lost or even reversed when the groups are combined.
Under H A , or if φ u = φ v , then we have: Proposition 2.1. Assume that the screening property holds and consider the sets I u , I v and I uv as fixed. Then under H A and regularity assumptions (A1)-(A3) (listed in Appendix A): A proof is given in Appendix A. Based on Proposition 2.1 we reject H 0 if LR nu,nv exceeds some critical value. The critical value is chosen to control the type-I error at some level of significance α. Alternatively, one may directly compute a p-value. In order to compute a critical value and/or determine a p-value we obtain in Section 2.3 the asymptotic distribution of LR nu,nv under H 0 .

Asymptotic Null Distribution
The work of Vuong (1989) specifies the asymptotic distribution of the log-likelihood ratio statistic for comparing two competing models in a very general setting, in particular Vuong's theory covers the case where the two models are non-nested. We apply Theorem 3.3 of Vuong (1989) to the special case where the competing models are M joint and M ind .
We further consider the matrices . The following theorem establishes the asymptotic distribution of LR nu,nv . Theorem 2.1 (Asymptotic null-distribution of restricted likelihood ratio-statistic.). Assume that the screening property holds and consider the sets I u , I v and I uv as fixed. If φ u = φ v (=:φ) and under regularity assumptions (A1)-(A6) (listed in Appendix A) we have: with r = |I u | + |I v | + |I uv | and Ψ r (·; ν) denotes the distribution function of a weighted sum of r independent chi-square distributions where the weights ν j , j = 1, . . . , r are eigenvalues of the matrix If we assume that the screening property holds then model M ind is correctly specified and the pseudo-true values (φ * u , φ * v ) equal the true values (φ u , φ v ). Furthermore, if we have φ u = φ v (=:φ) then the screening property guarantees that also M joint is correctly specified and that φ * uv =φ. We then write In Appendix A we prove the following proposition which characterizes the weights ν j defined as eigenvalues of the matrix W in Theorem 2.1. Proposition 2.2 (Characterization of eigenvalues). The eigenvalues ν j , j = 1, . . . |I u | + |I v | + |I uv |, of matrix W in Theorem 2.1 can be characterized as follows: • 2 |J| eigenvalues are 0.
• The remaining eigenvalues equal ± 1 − µ j , where µ j are eigenvalues of (2.8) Expressions (2.7) and (2.8) of Proposition 2.2 have nice interpretations in terms of analyzing the variances of the models M ind and M joint . Consider the random variables which are asymptotically Normal distributed with mean zero. If res c Ia and res c I b denote the residuals obtained from regressing Z c Ia and Z c I b against Z c J , then QI aIb is the covariance between these residuals. It is easy to see that the matrix QI uvIu Q − We further point out two special cases of Proposition 2.2: If M joint is nested in M ind , i.e., |I uv | = 0, then LR nu,nv follows asymptotically a chi-squared distribution with degrees of freedom equal |I u |+|I v |−|I uv |. If |I u | = |I v | = 0 then LR nu,nv is asymptotically distributed according to χ 2 |J| − χ 2 |Iuv|−|J| .

Estimation of the weights ν in Ψ r (·; ν)
In practice the weights ν of the weighted sum of chi-squared null distribution have to be estimated from the data. In light of Theorem 2.1, it would be straightforward to estimate the quantities B c IaI b (φ * a , φ * b ) and plug them into expression (2.6) to obtain an estimateŴ . Estimating ν then involves computation of r = |I u |+|I v |+|I uv | eigenvalues. Despite model reduction in the screening step, r can be a rather large number and can result in inefficient and inaccurate estimation. However, if both populations arise from the same distribution, then the overlap of the active-sets J = I uv ∩ I u ∩ I v is large compared to I u , I v and I uv . According to Proposition 2.2, the number of eigenvalues s = min{|I uv |−|J|, |I u |+|I v |−|J|} which remain to be estimated is small compared to r. We therefore estimate matrix (2.7) by and similarly matrix (2.8) by . One possibility is to use the sample analogues: Another way is to plug-in estimatorsφ u andφ v into the expectation with respect to φ u and φ v : (2.11) Figure 1 shows for a linear regression example with l = 100 predictors the estimated weightsν (upper panels showν obtained using directly Theorem 2.1; lower panels shoŵ ν obtained via Proposition 2.2). Estimating the eigenvalues with help of Proposition 2.2 works better than direct computation according to Theorem 2.1. In particular the true zero eigenvalues are poorly estimated with the direct approach. Figure 2 illustrates for the same example the quality of approximation of the ordinary and restricted likelihood-ratio with their asymptotic counterparts. The weighted sum of chi-squares approximates well already for small sample sizes, whereas the χ 2 p comes close to the distribution function of the ordinary likelihood-ratio only for very large n's.

P-Values and Multi-Splitting
In this section we demonstrate how to obtain p-values for testing hypothesis (1.1) using the methodology developed in Sections 2.1-2.4. The basic workflow is to split the data  into two parts, do screening on one part and derive asymptotic p-values on the other part. A p-value computed based on a one-time single split depends heavily on the arbitrary choice of the split and amounts to a "p-value lottery" (Meinshausen et al., 2009): for finite samples and for some specific split, the screening or the sparsity property might be violated which then results in a erroneous p-value. An alternative to a single arbitrary sample split is to split the data repeatedly. Meinshausen et al. (2009) demonstrated in the context of variable selection in the high-dimensional regression model that such a multi-split approach gives improved and more reproducible results. We adopt this idea and divide the data repeatedly into k = 1, . . . , K different splits. Let (U k in , V k in ) and (U k out , V k out ) be the first and second half of split k. On the first half screening is performed by solving three times the ℓ 1 -regularized log-likelihood problem (2.3), individually for each population and for both populations pooled together. The tuning parameter λ is always chosen by cross-validation. This gives models M k joint and M k ind defined via active-sets . Then, the restricted log-likelihood ratio LR k is evaluated on the second half of split k and a one-sided p-value is computed according to is defined in Theorem 2.1. The weightsν are estimated from the second half of split k as described in Section 2.4. As in Meinshausen et al. (2009) we aggregate p-values P k , k = 1, . . . , K, obtained from all different splits using the formula: where q γ (·) is the empirical γ-quantile function. This procedure is summarized in Algorithm 1. We refer to the choice K = 1 as the single-split method and to the choice K > 1 as the multi-split method.
Algorithm 1 Single-and multi-split algorithm for high-dimensional two-sample testing Input number of splits K. Randomly split data into (U k in , U k out ) and (V k in , V k out ).

Examples
As mentioned in the beginning of Section 2 two examples of our approach are differential regression where the populations follow linear regression models and differential network where the populations are generated from Gaussian graphical models. In this Section we provide details on both examples.

Differential regression Consider a regression model
with Y ⊂ R (k = 1), X ⊂ R l (l = p − 1) and random error ǫ ∼ N (0, σ 2 ). With φ = (β, σ 2 ) ∈ R p the score function is given by Now, given data U = (Y u , X u ) and V = (Y v , X v ) we obtain p-values by following Algorithm 1 (see Section 2.5). For the regression model, ℓ 1 -penalized maximum likelihood estimation, used in the screening step, coincides with the Lasso (Tibshirani, 1996) and is implemented in the R-package glmnet (Friedman et al., 2010). With the help of formula (3.13) we can easily compute plug-in estimatesB c IaI b ,plug-in (see equation (2.11)) and then the weightsν of the asymptotic null distribution are obtained as outlined in Section 2.4. We use the algorithm of Davies (1980), implemented in the R-package CompQuadForm (Duchesne and de Micheaux, 2010), to compute the distribution function of Ψ r (·;ν).
Differential network Let Y ⊂ R k (l = 0) be Gaussian distributed with zero mean and covariance Σ, i.e., N (0, Σ). A Gaussian graphical model with undirected graph G is then defined by locations of zero entries in the inverse covariance matrix Ω = Σ −1 , i.e., (j, j ′ ) ∈G ⇔ Ω jj ′ = 0. Setting φ = ({−0.5Ω jj } k j=1 , {−Ω jj ′ } j>j ′ ) ∈ R p (p = (k + 1)k/2), the score function is given by By invoking formulas on fourth moments of a multivariate normal distribution (see Appendix B) we find: In the Gaussian graphical model case ℓ 1 -regularized maximum likelihood estimation is well-known under the name Graphical Lasso (or GLasso) and is implemented in the Rpackage glasso (Friedman et al., 2008). As in differential regression, given data U = Y u and V = Y v , plug-in estimatesB c IaI b ,plug-in are computed by formula (3.14), weightsν are obtained subsequently as described in Section 2.4 and for p-value calculation we follow Algorithm 1.
Under H 0 , the regression coefficientsβ = β u = β v have 5 non-zero elements (at random location) with values 1.
Under H A , the regression coefficients β u and β v have entries of values 1 at three common locations and have entries of value α 1 at two different locations. α 1 = 0.25, 0.5.
Setting 2 As setting 1 but with SNR = 5.
Setting 3 As setting 1 but as a predictor matrix X we use for both populations gene expression data from n = 594 ovarian carcinomas. The data is publicly available at The Cancer Genome Atlas (TCGA) data portal (http://www.cancergenome.nih.gov). We select the l = 7, 10, 25, 50, 100, 250 genes exhibiting the highest empirical variance among samples. We choose σ u and σ v such that SNR=10.
Under H 0 , the inverse covariance matrices Ω u and Ω v are equal and have k non-zero entries at random locations.
The values α 1 and α 2 control the strength of the alternative for the regression and network cases respectively: A larger value α 1 results in two regression models which differ more. On the other hand a larger α 2 signifies that the two networks have more edges in common and are therefore more similar. For each of the four settings we perform 500 simulation runs for the null-and the different alternative hypothesis. We then compute for each run a p-value with the single-and the multi-split method (50 random splits, γ min = 0.05). We compare the proposed methodology with p-values obtained using the asymptotic χ 2 p -distribution of the ordinary likelihood-ratio statistic. We call this latter approach ordinary likelihood-ratio test. Further we compare with a permutation test where we use the symmetric Kullback-Leibler distance between ℓ 1 -regularized estimates as a test statistic. This permutation test uses 100 random permutations and further details are given in Appendix C. For each setting we evaluate the fraction of wrongly rejected null hypothesis (false positive rate) and the fraction of correctly rejected hypothesis (true positive rate) at a significance level of 5%.
Results are shown in Figures 3-6. The false positive rate (FPR) of the ordinary likelihoodratio explodes in all settings already for medium numbers of covariates l, respectively k. We further see that the single-split method is not able to control the true positive rate (TPR) at the 5% level in settings 1-2. This is most probably due to the fact that the sample size is too small and consequently the asymptotic approximation of the null distribution is not accurate enough. The multi-split method is much more conservative and has lower FPR. This observation is in agreement with results shown in Meinshausen et al. (2009). As expected, the permutation test controls the false positive rate at the 5% level. Concerning the TPR (or power of the test) we find that single-splitting, multi-splitting and the permutation test perform well in most settings. Only with weak alternatives (setting 1 and 2 with α 1 = 0.25) the TPR decreases. Interestingly, the permutation test exhibits loss of power for large values of l, whereas single-and multi-split display higher power than the permutation test for large l but lower power for small l. For completeness Figures 9-12 in Appendix D provide more insights on the single-split algorithm. These figures report results for settings 1-4 with respect to the sparsity and screening properties.
For setting 1 with p = 100 we further examine the distribution of the p-values obtained using the single-split method and the ordinary likelihood-ratio test when varying the sample size (n :=)n u = n v = 200, 500, 1000, 5000. For every n we use 100 samples in the screening step of the single-split method and the rest of the samples are used for pvalue calculations. Results are shown in Figure 7. We see that p-values obtained from the single-split method are very well approximated by the uniform distribution function. Only in the scenario with n = 200 the approximation is not accurate enough which is also reflected in a too large number of false positives. ["ordinary lrt": classical likelihood ratio test; "single-split", "multi-split": see text; perm: permutation test as described in text].

Application to genomic data from cancer biology
We apply the multi-split method to real datasets from cancer biology. For differential regression we take data from Broad-Novartis Cancer Cell Line Encyclopaedia (CCLE) (http://www.broadinstitute.org/ccle/home). Barretina et al. (2012) use the data to predict anticancer drug sensitivity from genomic data. They describe that the histone deacetylase inhibitor panobinostat shows increased sensitivity in haematological cancers ordinary lrt single−split multi−split perm Figure 5: Simulation study, Setting 3 (differential regression, real predictor matrix): False positive rate (FPR) and true positive rate (TPR) at 5%-level shown as a function of number of predictors (l). α 1 indicates the strength of the alternative (the larger α 1 the stronger the alternative). ["ordinary lrt": classical likelihood ratio test; "single-split", "multi-split": see text; perm: permutation test as described in text]. ordinary lrt single−split multi−split perm Figure 6: Simulation study, Setting 4 (differential network): False positive rate (FPR) and true positive rate (TPR) at 5%-level shown as a function of number of covariates (k). α 2 indicates the strength of the alternative (the smaller α 1 the stronger the alternative).
["ordinary lrt": classical likelihood ratio test; "single-split", "multi-split": see text; perm: permutation test as described in text]. ["ordinary lrt": classical likelihood ratio test; "single-split": see text]. compared to solid cancers. As the response variable we take experimentally determined sensitivity to panobinostat and we use gene expressions of the l = 100 genes showing highest Pearson correlation over all samples with panobinostat as the predictor matrix. We compare the three cancer subtypes with the largest sample size: lung cancer (89 cell lines), skin cancer (71 cell lines) and cancer with haematopoietic and lymphoid-tissue origin (abbreviated with haem, 40 cell lines). For differential network we use data from The Cancer Genome Atlas (TCGA) data portal (http://www.cancergenome.nih.gov). We consider gene expression data of the k = 42 genes present in the cancer gene list of Hahn and Weinberg (2002) (downloadable at http://www.cbio.mskcc.org/CancerGenes) and compare n u = 155 lung squamous cell carcinomas against n v = 174 colon adenocarcinomas.
For both examples we compute p-values with the multi-split method (500 splits) and with the permutation test (500 permutations). In the differential regression example the sample sizes are very small and cross-validation in the screening step can lead to activesets which are too large. Therefore, we adapt the screening step by setting the smallest coefficients to zero whenever there are more than ⌈5n⌉ non-zero coefficients (n denotes the sample size involved in the screening step). We further carry out "back-testing" by dividing data randomly into populations U and V and then performing significance testing. All results are shown in Table 1. In the regression example the multi-split method gives considerably smaller p-values than the permutation test. A potential explanation could be that the permutation test exhibits small power in difficult scenarios (see settings 1 and 2 in Section 4.1 with α 1 = 0.25, i.e., right panel in Figures 3 and 4). In the differential network example both methods have p-values which are numerically indistinguishable from zero. For the permutation test this implies that the observed test statistic is always larger than those obtained from random permutations. All p-values obtained for back-testing equal one which is reassuring. Figure 8 shows histograms of all 500 p-values obtained by the multi-split method and illustrates both the sensitivity of p-values to single splits and also the information contained in the entire distribution of p-values obtained in the iterative, multi-split approach. For example in the comparison of skin against haem we get the whole range of p-values between zero and one. However, the distribution of these p-values is heavily skewed towards zero which is reflected in an aggregated p-value of 0.022.

Discussion and conclusions
We have presented a novel and very general approach for high-dimensional two-sample testing. We combined ideas including sample-splitting, non-nested hypothesis testing and p-value aggregation to propose a methodology that allows two-sample testing for a wide range of high-dimensional models. We treated in detail linear regression (differential regression) and Gaussian graphical models (differential network ) and validated their performance on simulated and real data.
Our methodology is supported by asymptotic theory. However, our results obtained in Proposition 2.1 and Theorem 2.1 do not reflect a truly high-dimensional setup as they consider the active sets selected in the screening step as fixed. An aim of our ongoing research is to generalize these results to the case where the size of the active sets can grow with a rate of smaller order than the sample size.
Whilst the focus of this paper was not on applications, we note that the methodology we propose should have immediate utility in biology. Differential network can be used to test a number of hypotheses of current scientific interest. For example, in cancer biology it is widely believed that genomic differences between cancer types and subtypes may in some cases be manifested also at the level of biological networks, including gene regulatory and protein signaling networks. However, we are not aware of existing methodology that allows such hypotheses to be tested statistically under the relevant conditions of moderate sample size and high dimensionality. Our approach can be used to test such hypotheses directly from high-throughput data, as illustrated in the example above. A further application of the differential networks formulation of our approach is in gene-set testing. Currently, gene set tests (e.g., Subramanian et al. (2005); Irizarry et al. (2009)) are not truly multivariate. Our approach could be directly applied to test differences in gene sets at the level of not only means but also covariances or networks.
(A5) If (φ :=)φ u = φ v , thenφ is an interior point of Φ u , Φ v and Φ uv . Further,φ is a regular point of B M ind and B M joint .
(A6) The information matrix equivalence holds, i.e., for all A ⊂ {1, . . . , p} we have: Proof of Proposition 2.1 Given assumptions (A1)-(A3) and noting that φ * is the standard Kullback-Leibler divergence. If φ u = φ v , then at least one term on the right hand side of equation (A.15) is strictly positive. Therefore, under H A : LR n,n a.s.

→ ∞ (n → ∞).
Proof of Theorem 2.1 Consider the setting in Vuong (1989) with competing models F φ = M ind and G γ = M joint , i.e., f (y|x; θ) = d(y u |x u ; φ u )d(y v |x v ; φ v ) and g(y|x; γ) = d(y u |x u ; φ uv )d(y v |x v ; φ uv ) with y = (y u , y v ), x = (x u , x v ), θ = (φ u , φ v ) and γ = φ uv . If the screening property holds then model M ind is correctly specified and the pseudo-true values (φ * u , φ * v ) equal the true values (φ u , φ v ). Furthermore, if we have (φ :=)φ u = φ v then the screening property guarantees that also M joint is correctly specified and that φ * uv =φ. As a consequence we have f (y|x; θ * ) = g(y|x; γ * ). Therefore, assuming assumptions (A1)-(A5), it follows from Theorem 3.3 (i) in Vuong (1989) that where r = |I u | + |I v | + |I uv | and the weights ν j , j = 1, . . . , r, are eigenvalues of the matrix The matrices B f , B g , B f g , B gf , A f and A g are defined as in Vuong (1989). If we set then we obtain as a consequence of the independence of the populations U and V : IuvIv . By invoking the information matrix equivalence, assumption (A6), we find Proof of Proposition 2.2 Set r u = |I u |, r v = |I v |, r uv = |I uv | and r = r u + r v + r uv . The eigenvalues of W are the solutions to the equation If ν = 1, then we find by using Schur's complement where µ = (1 − ν 2 ). The second term on the right of equation (A.17) has r uv roots. Therefore r u + r v − r uv of the eigenvalues ν j equal one. As φ * u = φ * v = φ * uv , we can write An easy calculation involving Schur's complement shows Further, we find Putting together, equation (A.18) reads and conclude that 2|J| of the eigenvalues are zero and that the remaining eigenvalues are obtained by computing eigenvalues of The proof for the case with |I uv | > |I u | + |I v | works similarly.

C Permutation test and symmetric Kullback-Leibler divergence
The symmetric Kullback Leibler divergence is defined as Consider a random partitition of the data into two groupsŨ = (Yũ, Xũ) andṼ = (Yṽ, Xṽ). We construct a permutation test based on the test statistic For the regression model the Kullback-Leibler divergence equals: For the Gaussian graphical model we have: Figure 10: Simulation study, Setting 2 (differential regression, SNR=5): Same notation as in Figure 9.  Figure 11: Simulation study, Setting 3 (differential regression, real predictor matrix): Same notation as in Figure 9.  Figure 12: Simulation study, Setting 4 (differential network): Same notation as in Figure 9.