{'fraction_non_alphanumeric': 0.047264954691874646, 'fraction_numerical': 0.014841797874582292, 'mean_word_length': 4.68120403625791, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiple multi-sample testing under arbitrary covariance dependency

Modern high-throughput biomedical devices routinely produce data on a large scale, and the analysis of high-dimensional datasets has become commonplace in biomedical studies. However, given thousands or tens of thousands of measured variables in these datasets, extracting meaningful features poses a challenge. In this article, we propose a procedure to evaluate the strength of the associations between a nominal (categorical) response variable and multiple features simultaneously. Specifically, we propose a framework of large-scale multiple testing under arbitrary correlation dependency among test statistics. First, marginal multinomial regressions are performed for each feature individually. Second, we use an approach of multiple marginal models for each baseline-category pair to establish asymptotic joint normality of the stacked vector of the marginal multinomial regression coefficients. Third, we estimate the (limiting) covariance matrix between the estimated coefficients from all marginal models. Finally, our approach approximates the realized false discovery proportion of a thresholding procedure for the marginal p-values, for each baseline-category pair. The proposed approach offers a sensible trade-off between the expected numbers of true and false rejections. Furthermore, we demonstrate a practical application of the method on hyperspectral imaging data. This dataset is obtained by a matrix-assisted laser desorption/ionization (MALDI) instrument. MALDI demonstrates tremendous potential for clinical diagnosis, particularly for cancer research. In our application, the nominal response categories represent cancer subtypes.


INTRODUCTION
Many datasets in various scientific disciplines, such as neuroimaging, genomics, brain-computer interfacing, and others, are nowadays high-dimensional. Such datasets are frequently analysed by means of large-scale multiple testing for some phenotype -for example, cancer subtypes -on each feature among thousands of features. This approach has multiple applications in the life sciences. 1 At least since the seminal publication 2 , control of the false discovery rate (FDR) is a commonly used type I error criterion in high-dimensional multiple test problems. However, applying the well-known Benjamini-Hochberg 2 (henceforth, B-H) procedure or Storey's 3,4 procedure, which have originally been designed for stochastically independent or weakly dependent p-values, respectively, can lead to an FDR inflation under certain forms of (strong) dependencies. 5,6,7,8 Under the assumption of positive regression dependency on subsets (PRDS), the B-H procedure still controls the FDR 9 , but checking the validity of the PRDS assumption on the basis of a sample is far from trivial and only established for a limited number of special cases like, for example, elliptically distributed vectors of test statistics. 10 For these reasons and with the goal of optimizing the statistical power of the multiple test under the constraint of (at least asymptotic) FDR control at the desired level, multivariate (FDR-controlling) multiple tests have drawn a lot of attention over recent years. 11 A multivariate multiple test incorporates the dependence structure among test statistics or p-values, respectively, or an estimate thereof, explicitly in its decision rule. Large-scale multiple testing under general and strong dependency remains challenging and an active research topic in modern statistics. 7,8,12,13 Among other approaches like, for instance, considering block dependencies in genetics 14,15 , several multi-factor models have been proposed to model dependencies among test statistics. 7,8,13,5,6 In particular, a general framework to approximate the false discovery proportion (FDP) has been introduced. 7, 16 It employs a so-called principal factor approximation (PFA) of the (estimated limiting) covariance matrix of test statistics, assuming that these test statistics are (at least asymptotically for large sample sizes) jointly normally distributed. The general strategy of the approach is to apply a spectral decomposition of the aforementioned covariance matrix and then to deduct its principal factors that mostly induce the strong correlation dependency. Since the FDR is the expected value of the FDP, bounding the FDP (with high probability) ensures FDR control, too.
Biological datasets generated by advanced high-throughput devices typically contain thousands of measured variables, many of which are related. In this work, we consider a dataset obtained by a matrix-assisted laser desorption/ionization (MALDI) imaging mass spectrometry (IMS) tool, also known as MALDI imaging. MALDI has advanced considerably and demonstrates immense potential in numerous pathological applications. 17,18 MALDI-generated datasets contain thousands of measured variables (corresponding to molecular masses), many of which are (highly) related. 19 Since MALDI data do in general not exhibit dependencies in blocks, it has recently been advocated to model dependencies among MALDI features by means of factor models, and to employ PFA in the context of multiple two-sample comparisons based on MALDI data. 20 In this work, we deal with multiple multi-sample hypotheses testing under arbitrary correlation dependency, in particular in the case of more than two samples or categories, respectively. We do so under the scope of multiple marginal multinomial regression models and by applying the PFA methodology. We call our method multinomial -principal factor approximation ("Multi-PFA").
The remainder of this article is structured as follows. In Section 2, we present the proposed multi-sample multiple testing approach. Section 3 is dedicated to a simulation study under different data-generating scenarios. We demonstrate the practical application of the proposed approach to hyperspectral data in Section 4, and we conclude with a discussion and an outlook in Section 5.

Data Structures
MALDI imaging data are typically stored in an × matrix , where rows correspond to mass spectra and columns correspond to thousands of mass-to-charge (m/z) values. 21 Since MALDI data represent molecular masses of ionizable molecules, ∈ ℝ × ≥0 . However, our proposed framework makes no explicit assumption that data points have to be non-negative. We approach the biological challenge of analyzing the associations between (individual) m/z values (predictors ( ) 1≤ ≤ ) and several cancer subtypes (categorical response ) by carrying out multiple statistical hypothesis tests simultaneously. Under the -th null hypothesis 0 , there is no association between and , and we aim at testing the null hypotheses 01 , … , 0 (against their two-sided alternatives 11 , … , 1 ) simultaneously based on one and the same dataset. Throughout the remainder, we assume that is a (random) nominal outcome variable, meaning that the response categories are unordered. Hence, for each ∈ {1, … , }, the tuple ( , ) takes its values in ℝ ≥0 ×{1, 2, … , }, where we denote by the number of (response) categories.

Marginal Modelling for Categorical Responses
First, we model marginal associations for each predictor individually. The reason for considering marginal modelling relies on the fact that we consider an inferential framework, in which the number of covariates is either quite large or even (much) larger than the number of observational units ( << ). Thus, a model that includes all predictors in a single model cannot be fitted (reliably). In the context of large-scale multiple (association) testing, marginal linear regression models have been adopted to test an individual hypothesis for each predictor. 22,7 Marginal generalized linear models (GLMs), for positive outcome variables, have been considered in situations where the focus of the statistical inference is on the mean of the response given a group of predictors. 23 Also, marginal logistic regressions have been utilized to analyse child obesity 24 , with the usage of generalized estimating equations for jointly estimating (univariate) associations. 25 For this section, let ∈ {1, … , } be arbitrary, but fixed. To model the conditional distribution ℙ( | ), we assume a marginal multinomial regression with a canonical (logit) link function. In general, the response categories cannot be treated as a one-dimensional response, and we have to introduce a dummy variable for each but one category. Therefore, we have a multivariate outcome variable in the (relevant) case of > 2 (for more details, see Chapter 7 in 26 ). One (response) category is used as a baseline, and the multinomial logistic model is established by pairing each of the other nominal response categories with the baseline. Throughout this article, response category is used as the baseline. Furthermore, we mainly consider the case of = 3 throughout the remainder, for concreteness and because of its relevance for the application that we are going to present in Section 4. For each (marginal) fit, we consider intercept terms and slope parameters , for 1 ≤ ≤ − 1. We denote the observables (regressor and response variable) for the -th marginal regression fit by ( ( ) , ( ) ) 1≤ ≤ , and we assume that these tuples are stochastically independent and identically distributed bivariate random vectors, each having the same (joint) distribution as ( , ). For observational unit ∈ {1, … , }, our model equation for the -th baseline-category logit is given by Substituting log( ( ) ∕ ( ) ) = + ( ) as well as ( ) = 1∕{1 + ∑ −1 =1 exp( + ( ) )} in (2), we find that where stands for the log-likelihood function, and̂ and̂ denote the ML estimators. In practice, the estimation is carried out conditionally to the actual observed values of the response indicators and the regressors, leading to the numerical values of the ML estimates.
In this work, we are interested in testing simultaneously two families of hypothesis-alternative pairs, which are given by for = 3. By means of testing the hypotheses in (4), we make binary decisions (rejection or non-rejection) for both baselinecategory sets, consisting of null hypotheses each, based on the data at hand. In the MALDI context, this means that we aim at identifying the most distinctive m/z values for the cancer association, where three different cancer types are considered in our application.

Multiple Marginal Models
The following (second) step of the proposed framework is to combine all marginal fits and to approximate the joint null distribution of all ML estimators within each baseline-category pair. To do so, we exploit the approach for jointly estimating the parameters of multiple marginal models developed by Pipper et al. 27 , and apply this approach to the marginal models explained in the previous section. Notice that our modelling approach implies that the regression parameters are unique to one model and are not shared by any two models 1 ≠ 2 . In addition to that, the intercepts ( ) 1≤ ≤ , 1≤ ≤2 are nuisance parameters in the sense that the hypotheses in (4) solely concern the 's. However, the intercepts ( ) 1≤ ≤ , ∈{1,2} contribute to the estimation and the standardization of 's. We aim at establishing a central limit theorem for the two -dimensional random vectorŝ = (̂ 1 , … ,̂ ) ⊤ 1≤ ≤2 . This is accomplished by stacking the standardised score contributions of thê 's across all marginal models, within each baselinecategory pair. To this end, we use that, under standard regularity assumptions (like finite variances and non-vanishing limiting relative category frequencies), each ML estimator̂ admits the asymptotic (i. e., → ∞) representation 27 is the relevant row of the inverse Fisher information matrix of the model specified in (1) that corresponds to coordinate for the -th observational unit,Ψ is the score function pertaining to coordinate and category for the -th observational unit, and ℙ (1) indicates a sequence of random variables converging to zero in probability.
Elementary calculations yield that (Ψ 1 , Ψ 2 ) ⊤ is given by the second and the fourth coordinate of the four-variate vector where and ⊗ denotes the Kronecker product. Moreover, the Fisher information matrix of the full sample of size is the sum of the Fisher information matrices pertaining to each ∈ {1, … , }, due to the assumption of stochastically independent observational units. For more (computational) details regarding multinomial regression, see, for example, 28,29 .
and ∈ {1, … , }. Then, we can conclude from (5) the asymptotic expansion where now ℙ (1) indicates a sequence of -dimensional random vectors converging to the zero vector in probability. Since the observational units are assumed to be stochastically independent and identically distributed, the multivariate central limit theorem yields convergence in distribution (indicated by the symbol ← ← ← ← ← ← →) to a centered -variate normal distribution, i. e., that For each ∈ {1, 2}, the limiting variance-covariance matrix Σ can be estimated consistently bŷ where eachΨ is obtained by inserting the parameter estimates (̂ ,̂ ) 1≤ ≤ , ∈{1,2} from the marginal fits (cf. (3)) instead of the true parameter values ( , be the standardised versions of̂ 1 , …̂ , which are given by where √V ar(̂ ) is the square root of the -th diagonal element ofΣ , divided by √ . We have that pertaining toΣ from (9), and the notation ∼ approx.
indicates the approximate distribution for large . Therefore, since √V ar(̂ ) > 0 holds true with probability one, the families of hypotheses from (4) can be equivalently written as For the sake of presentation, we use henceforthΣ * and 1 , … , only for the first baseline-category pair (i. e., for = 1) and omit the index , for convenience of notation. However, we consider the same inference procedure for both baseline-category pairs.

False Discovery Proportion
Considering the multiple test problem given by (12) (for = 1), let 0 = |{ ∶ = 0}| denote the number of true null hypotheses, and 1 = − 0 = |{ ∶ ≠ 0}| the number of false null hypotheses. For controlling type I errors, we employ empirical process techniques, in the spirit of Storey's method. 3 In particular, we consider thresholding rules which reject a null hypothesis 0 if and only if a corresponding -value is smaller than or equal to a (data-dependent) threshold . Such types of multiple tests have been considered widely in previous literature. 7,16,3 Conceptually, the goal of the proposed method is to approximate the proportion of false discoveries among all rejections (commonly referred to as the false discovery proportion, FDP for short) for a fixed threshold under an arbitrary correlation matrix Σ * . To this end, we consider the three empirical processes where ranges in [0, 1] and the notation is used to indicate that we consider the -th -value as a random variable here. For fixed ∈ [0, 1], the quantities ( ), ( ), and ( ) represent the (random) number of false discoveries (i. e., false rejections or, synonymously, type I errors), the (random) number of true discoveries, and the (random) total number of discoveries, respectively, hence ( ) = ( ) + ( ). The FDP is for a fixed given by FDP( ) = ( )∕ max{ ( ), 1}, where the maximum is taken to avoid an expression of the form 0∕0. All of the latter random variables depend on the -statistics 1 , 2 , … , , because each (random) -value will be obtained by a transformation of , for 1 ≤ ≤ , see below. Notice also that ( ) and ( ) are both unobservable, whereas ( ) is observable.

Principal Factor Approximation
The next step of our proposed framework is to incorporate the correlation effects of the -statistics in an approximation of FDP( ) for a fixed threshold . The procedure relies on the identification of a low-dimensional linear space of random vectors which captures most of the dependence structure of the -statistics.
To this end, we carry out the spectral decomposition of the estimated correlation matrixΣ * . We let 1 , … , ≥ 0 denote the eigenvalues ofΣ * , and 1 , … , the corresponding eigenvectors. Hence,Σ * is represented by its eigenvalue-eigenvector pairs ( , ) 1≤ ≤ . For a fixed number of common factors, where 1 ≤ << , the factor approximation ofΣ * is given bŷ ). If is chosen appropriately, the × -matrix , corresponding to latent variables, can describe most of the dependence structure among the -statistics. Accordingly, 1 , … , can be decomposed as where is the vector of (latent) common factors. Due to orthogonality projection, these common factors are jointly stochastically independent. The random vector ( 1 , … , ) ⊤ ∼ (0, ) is called the vector of random errors, and it is assumed that factors and random errors are stochastically independent. In our previous work 20 , we elaborated more on the importance and the choice of .
On the basis of the above considerations, we consider the FDP estimator from Proposition 2 by Fan et al. 7 , which is for a given ∈ [0, 1] defined asF In this, the (unadjusted) two-sided, random -value corresponding to is given by = 2Φ(−| |) = 2(1 − Φ(| |)), and = ∑ ℎ=1 ℎ̂ is a linear estimator for = ⊤ . The FDP estimator given in (15) is based on a sparsity assumption. 7,16 Specifically, it is assumed that the numbers and 0 are large, while the number 1 of false nulls is relatively small. This assumption allows for the summation over all hypotheses (not only over "true nulls") in the numerator of (15). In the context of MALDI modelling, this sparsity assumption is justified, for three reasons: First, several researchers who have analysed MALDI datasets have reported a low number of m/z values that are highly associative for a cancer subtype. 30,18 Second, to the best of our knowledge, the reported number of biologically meaningful molecules, so-called biomarkers, in MALDI-related studies, is low -from three to five per study. 31,19 Third, there is evidence of a small signal-to-noise ratio 19,18 .
To evaluate (15) in a practical application, one needs to determine the (linear) estimator̂ = (̂ 1 , … ,̂ ) ⊤ of the common factors. The authors of 7 have proposed to estimatê via a linear regression ( 2 -regression) or by using a quantile regression ( 1 -regression). Regarding the 2 estimator, it has been proposed to carry out the estimation of (̂ 1 , … ,̂ ) ⊤ only on the basis of the subset of length 0.95 of the smallest (in absolute values) -statistics. This leads to the 2 estimator where the Z-values in (16) are sorted in ascending order based on their absolute values. We have utilized this estimator in our simulation study (see Section 3). The estimator based on the 1 -regression is given bŷ For the analysis of the MALDI data (see Section 4), we have used the 1 -regression rather than the 2 -regression, because it is more robust to highly untypical observations (i. e., outliers). Finally, the dependency-adjusted (random) -values pertaining to 1 , … , are given bỹ For a given threshold ∈ [0, 1], the null hypothesis 0 from (12) gets rejected based on the data at hand, if and only if̃ ≤ , 1 ≤ ≤ . For the purpose of FDP control, can be selected as the largest value = ∈ [0, 1] fulfilling thatFDP( ) is not exceeding a pre-determined level ∈ (0, 1), for example, = 10%. In practice, one can carry out a grid search over a set of candidate values of in order to find the value for a given . 20

Simulation Setup
We carried out a simulation study to assess the performance of our proposed methodology under different data-generating schemes. We considered the parameter settings = 500, ∈ {500, 1000}, and 1 = 10 (for each baseline-category pair). For each combination of the aforementioned parameter values, we performed 1, 000 Monte Carlo repetitions. Without loss of generality, we set 1 ≠ 0 as well as 2 ≠ 0 for ∈ {1, … , 1 }, and we call these first 1 coordinates "active". The remaining 0 coordinates have been set "inactive", meaning that = 0 for all in ∈ { 1 + 1, … , } and ∈ {1, 2}. Furthermore, we set all intercept terms ( ) 1≤ ≤ , ∈{1,2} to zero in all simulations. We estimated the common factors { ℎ ∶ ℎ = 1, … , } by applying the least-squares estimator as given in (16). The utilized number of common factors is reported in each caption of Tables 1 -4 below.
For each observational unit ∈ {1, … , }, the simulation data have been generated according to the model The considered dependency structures for 1 , … , are given in Model 1.
Scenario 2: 1 +1 , … , are jointly normally distributed on ℝ 0 , namely, according to 0 (0, Σ * ), where Σ * is an equicorrelation matrix with diagonal elements equal to one, and off-diagonal elements equal to some given value . Furthermore, the subvector ( 1 +1 , … , ) ⊤ is stochastically independent of the subvector ( 1 , … , 1 ) ⊤ . This is in order to avoid spurious effects of covariates with 1 + 1 ≤ ≤ on the response, which would arise from confounding if the two aforementioned subvectors would be dependent. The random variables 1 , … , 1 are stochastically independent and identically (0, 1)-distributed. The equi-correlation model is a one-factor model with a dominating first eigenvalue; see, e. g., Example 2.1 in 32 .
This simulation study has been run in R version 4.1.1, 33 employing the R function rmvnorm() 34 to simulate correlation dependency among the independent variables. In addition, we have employed the function vglm() 35 for the estimation of the marginal multinomial estimates. The function pfa.test() 36 has been used for applying the PFA method.

Simulation Results
We summarize the results from our computer simulations in Tables 1 -4 in terms ofFDP( ), ( ), and ( ) for a fixed threshold = 10 −4 . In addition, we report the median value of for the conventional choice of = 0.05 as well as the average (over the 1,000 Monte Carlo repetitions) of ( ). Due to the choice of 1 = 10, the largest possible value of ( ) equals ten, for any choice of the threshold .
Tables 1 and 2 contain our simulation results under Scenario 1 from Model 1. Because the -statistics are jointly stochastically independent (due to joint independence among all ), 0.05 is rather small here, because the "effective number of tests" (in the sense of Section 3.4 in 11 ) equals under joint independence of all test statistics or -values, respectively. This means that a rather strong multiplicity correction is necessary. However, the standard error ofFDP( ) is also quite small under Scenario 1, since the FDP is well concentrated around its expectation (the FDR) under joint independence among all test statistics orvalues, respectively. With regard to type II errors, the reported average values of ( 0.05 ) demonstrate that on average 7 − 8 out of 10 active coordinates could be identified by our proposed multiple test procedure.
Simulation results under Scenario 2 (Gaussian equi-correlation model) are summarised in Tables 3 and 4. Here, the effective number of tests is smaller than whenever ≥ 0, and it decreases as increases. Consequently, 0.05 grows with , too. In turn, this also leads to an improved power of the proposed multiple test, which is reflected by the numbers reported for ( ), which are under Scenario 2 on average larger than the corresponding values under Scenario 1.

Description of the Dataset
We applied our proposed inferential procedure to a hyperspectral dataset obtained from a MALDI imaging instrument. For a detailed description of this dataset in terms of its acquisition protocols, tissue sections, tissue blocks, etc., see 30,18 . In the aforementioned studies, the researchers analysed this dataset by merging two lung cancer subtypes into a single category. By doing so, the authors worked with two cancer subtypes based on two distinctive human organs (the so-called "LP task", i.e. lung vs. pancreas). In our study, however, we model this dataset with a three-class model, using the subtype of pancreatic The total number of hypotheses equals = 500; the number of false null hypotheses equals 1 = 10 per baseline-category pair (non-zero regression coefficients equal 1); the number of factors equals = 10; the rejection threshold equals = 10 −4 , apart from the last two columns. The total number of hypotheses equals = 1000; the number of false null hypotheses equals 1 = 10 per baseline-category pair (non-zero regression coefficients equal 1); the number of factors equals = 10; the rejection threshold equals = 10 −4 , apart from the last two columns.
adenocarcinoma as the baseline. The resulting goal of the statistical analysis is to perform simultaneously two (sub-) tasks, namely the pancreatic vs. lung adenocarcinoma and the pancreactic vs. lung squamous cell carcinoma comparison. As mentioned in the introduction, the high-throughput device that motivated this study is MALDI IMS. In general, this technology provides molecular information about a given analyte (for example, tissue) in a spatial manner. More concretely, the MALDI IMS tool measures mass spectra at multiple discrete spatial positions and yields an image for each unique spot within a provided tissue. For an illustration, Figure 1 displays instances of outputs from a MALDI experiment from three spatial spots. In Figure 1, the m/z values are plotted on the horizontal axes, and the relative abundances (intensities values) of ionizable molecules are plotted on the vertical axes. This "spatial molecular (biochemical) information can then be used for the determination of the cancer subtypes or the identification of the origin of the primary tumour in patients" 30 . Each spot is called a mass spectrum 37,21 , and it depicts relative abundances of ionizable molecules with numerous mass-to-charge (m/z) values -"ranging from several hundred up to a few tens of thousands m/z values" 21 -whereas a single m/z value in the context of MALDI is interpreted as a molecular mass. The full pipeline of a MALDI experiment, from a tissue spot to data analysis, is thoroughly illustrated in Figure 1 of 18 and Figure C2 of 19 . In the dataset analysed here, formalin-fixed paraffin-embedded (FFPE) 38 tumour samples from 445 patients were provided by the tissue bank of the National Center for Tumor Diseseas (NCT) Heidelberg, Germany. Tissue cores of three cancer subtypes were used, comprising lung adenocarcinoma (ADC) (168 patients), lung squamous cell carcinoma (SqCC) (136 patients) and pancreatic adenocarnimoma (141 patients). The corresponding dataset is publicly avaible from (https://gitlab.informatik.uni-bremen.de/digipath/Deep_Learning_for_Tumor_Classification_in_IMSlast accessed 7 February 2022).

Data Preprocessing
We performed the following (further) data preprocessing steps. First, due to the high dimensionality of the feature space (27, 286 m/z values), spectral filtering was carried out. Namely, we loaded the data into MATLAB R2018, and used the internal library The total number of hypotheses equals = 500; the number of false null hypotheses equals 1 = 10 per baseline-category pair (non-zero regression coefficients equal 1); the number of factors equals = 1; the rejection threshold equals = 10 −4 , apart from the last two columns. The total number of hypotheses equals = 1000; the number of false null hypotheses equals 1 = 10 per baseline-category pair (non-zero regression coefficients equal 1); the number of factors equals = 1; the rejection threshold equals = 10 −4 , apart from the last two columns.
MSClassifyLib (function MSAdaptiveResampleTrigger()) to apply the spectral filtering with default value of 0.4 Da (for more details, see, e.g., 19,18 and references therein). Second, the mass range was pruned up to the mass range of 2100 m/z, and we used one mass spectrum from each patient. Finally, the well-established normalisation step TIC (total ion count) 37 was employed on the chosen mass spectra range. The result of this data preprocessing was a dataset with 445 mass spectra (observational units) and 1579 m/z channels (features).

Analysis of the MALDI Dataset
We have modelled each pixel (based on the averaging filter around 0.4 Da as explained in the previous section) individually to determine those m/z values that are highly associative for a certain cancer subtype. Thus, the m/z values take the role of the 's and the cancer subtype takes the role of from our general setup, where ∈ {1, … , = 1,579} here. Statistically speaking, our procedure identifies significant differences in distributions across the subtypes, and is not only identifying significant peaks in the spectra. Figure 2 displays the empirical distributions of the Z-statistics, along with histograms of the non-adjusted -values, for both tasks. Clearly, the -values (for both tasks) do not resemble the realization of a random sample from the theoretical null distribution, which is the standard normal distribution on ℝ. In particular, the empirical variances are much larger than one in both histograms displayed in Figure 2. There are two reasons for this overdispersion, namely, (i) the presence of correlation effects among the -values (for more details, see 5,6 ) and (ii) the presence of extremely large (in absolute value) -values, which is likely caused by significant effects, peaks and isotopic patterns. In our previous study 20 , we have observed a similar phenomenon and have described a procedure (based on the "empirical null distribution") which can be used to account for the part of the overdispersion which is caused by correlations among the -statistics.
The main results of our real data analysis are presented in Figure 3, where we have chosen = 3 for both tasks. According to the number = 1,579 and taking into account the magnitude of the observed correlation effects among the -statistics, a plausible range for the rejection threshold has been selected for being displayed on the horizontal axes. The subplots in Figure 3 illustrate the total number of rejections ( ( )), the estimated number of false discoveries ( ( )) and the estimated FDP (FDP( )) over this range. Apparently, all three aforementioned quantities decrease with decreasing (notice the negative logarithmic scale of the horizontal axes in Figure 3). For both tasks,FDP( ) ranges in [6%, 17%] for ∈ [10 −7 , 10 −5 ]. Table 5 tabulates ( ) and FDP( ) for several threshold values from the latter range.
To the best of our knowledge, no biomarkers have been confirmed yet for this dataset (pancreas vs. lung task). 30 However, according to the previous reports 30,18 Tables 6 and 7 indicate those m/z values that are distinctive for the pancreatic cancer subtype.

DISCUSSION AND OUTLOOK
Motivated by MALDI association studies, this study's objective has been to evaluate the strength of associations between a nominal variable of interest -describing cancer subtypes in our application -and a large number of measured variables (in our case given by m/z values). We have proposed an approach to screen all features and identify the most associative ones with the multi-class outcome. Our approach decomposes the outcome variable into multiple baseline-category pairs, and approximates  the false discovery proportion under arbitrary correlation dependency within each pair. As demonstrated on simulated data, this approach leads to a sensible balance between the (average) numbers of true and false rejections. Furthermore, both our simulation study and the presented application to real data demonstrate that the proposed method can be applied even in cases where the sample size is considerably smaller than the number of features. This makes the procedure attractive for medical screening applications. Two central assumptions underlie the proposed methodology: First, we assume sparsity in the sense that the number of active features is small. Second, we assume that the dependency structure amongst the considered -statistics can accurately be approximated by a multi-factor model. The approach of MMM, however, is flexible in the sense that it does not depend on heavy assumptions. For its applicability, we only have to assume that the multi-class outcome is associated with features (m/z values), and that this association can accurately be described by a (marginal) multinomial regression model for each m/z-value separately.  From the application perspective, we have applied the Multi-PFA method to a MALDI imaging dataset consisting of a large number of m/z values and one spectrum from each patient. To the best of our knowledge, we are the first to address the threeclass problem for this dataset. In this way, we have provided a more detailed analysis for both lung cancer subtypes than previous studies.
There are several potential directions for further research. First, it might be captivating to consider different supervised learning methods (for instance, a neural network with more than one layer) instead of the multinomial regression model proposed in this paper. Second, it may be of interest to evaluate the uncertainty regarding the realized FDP in order to provide a confidence region for it, in addition to mere point estimation of the FDP. Third, disentangling correlation effects and regression effects on the empirical -statistics distribution in a more detailed manner is an interesting follow-up research topic in our context.