Meaningful Problem Instances and Generalizable Results

Schriftenreihe CIplus, Band 1/2015 Herausgeber: T. Bartz-Beielstein, W. Konen, H. Stenzel, B. Naujoks Meaningful Problem Instances and Generalizable Results Thomas Bartz-Beielstein Meaningful Problem Instances and Generalizable Results∗ Thomas Bartz-Beielstein SPOTSeven Lab Cologne University of Applied Sciences Steinm¨ ullerallee 1 51643 Gummersbach [email protected] www.spotseven.de February 12, 2015 Abstract Computational intelligence methods have gained importance in several real-world domains such as process optimization, system identification, data mining, or statistical quality control. Tools are missing, which determine the applicability of computational intelligence methods in these application domains in an objective manner. Statistics provide methods for comparing algorithms on certain data sets. In the past, several test suites were presented and considered as state of the art. However, there are several drawbacks of these test suites, namely: (i) problem instances are somehow artificial and have no direct link to real-world settings; (ii) since there is a fixed number of test instances, algorithms can be fitted or tuned to this specific and very limited set of test functions; (iii) statistical tools for comparisons of several algorithms on several test problem instances are relatively complex and not easily to analyze. We propose a methodology to overcome these difficulties. It is based on standard ideas from statistics: analysis of variance and its extension to mixed models. This chapter combines essential ideas from two approaches: problem generation and statistical analysis of computer experiments. ∗This is a preprint of the publication T. Bartz-Beielstein. How to create generalizable results. In J. Kacprzyk and W. Pedrycz, editors, Springer Handbook of Computational Intelligence, chapter 56. Springer, 2015 (in print). The original publication is available at www.springerlink.com 1 1 Introduction Computational intelligence (CI) methods have gained importance in several realworld domains such as process optimization, system identification, data mining, or statistical quality control. Tools are missing, which determine the applicability of CI methods in these application domains in an objective manner. Statistics provide methods for comparing algorithms on certain data sets. In the past, several test suites were presented and considered as state of the art. However, there are several drawbacks of these test suites, namely: •problem instances are mostly artificial and have no direct link to real-world settings; •since there is a fixed number of test instances, algorithms can be fitted or tuned to this specific and very limited set of test functions. As a consequence, studies (benchmarks) provide insight how these algorithms perform on this specific set of test instances, but no insight on how they perform in general; •statistical tools for comparisons of several algorithms on several test problem instances are relatively complex and not easy to analyze. We propose a methodology to overcome these difficulties. This methodology, which generates problem classes rather than uses one instance, is constructed as follows. 1. First, we pre-process the underlying real-world data. 2. In a second step, features from these data are extracted. This extraction relies on the assumption that mathematical variables can be used to represent real-world features. For example, decomposition techniques can be applied to model the underlying data structures, if we are using time-series data. The original time series is deconstructed into a number of component series, where each of these reflects a certain type of behavior, e.g., a trend or seasonality[9]. We obtain an analytic model of the data. 3. Then, we parametrize this model. Based on this parametrization and randomization, we can generate infinitely many new problem instances. 4. If no real-world data are available, problem instances can be generated using test-problem generators. The generation of test problems, which are well-founded and have practical relevance, is an on-going field of research for several decades. 5. From this infinite set, we can draw a limited number of problem instances which will be used for the comparison. 6. Since problem instances are selected randomly, we apply random and mixed models for the analysis [15]. Mixed models include fixed and random effects. A fixed effect is an unknown constant. Its estimation from the data is a common practice in analysis of variance (ANOVA) or regression. A random effect is a random variable. We are estimating the 2 parameters that describe its distribution, because—in contrast to fixed effects—it makes no sense to estimate the random effect itself. This chapter combines ideas from two approaches: problem generation and statistical analysis of computer experiments. The work presented by Chiarandini and Goegebeur [11] provides the basis of our statistical analysis. They present a systematic and well-developed framework for mixed models. Related modeling approaches were suggested by McGeoch[14] and Birattari [7]. Gallagher and Yuan [13] present a problem instance (landscape) generator that is parameterized by a small number of parameters, and the values of these parameters have a direct and intuitive interpretation in terms of the geometric features of the landscapes that they produce. Casti˜neiras, Cauwer, and O’Sullivan [10] present a parametrizable benchmark generator for bin packing instances based on the well-known Weibull distribution. Using the shape and scale parameters of the Weibull distribution, the authors generate benchmarks that contain a variety of item size distributions. They report that for all bin capacities, the number of bins required in an optimal solution increases as the Weibull shape parameter increases. Using this feature, scalability is enabled. Basically, this chapter tries to find answers for the following fundamental questions in experimental research. (Q-1) How to generate problem instances? (Q-2) How to generalize experimental results? The chapter is structured as follows. Section 2 introduces real-world and artificial optimization problems. Algorithms are described in Sect. 3. Objective functions and statistical models are introduced in Sect. 4. These models take problem and algorithm features into consideration. Section 5 presents case studies, which illustrate our methodology. This chapter closes with a summary and an outlook. 2 Features of Optimization Problems 2.1 Problem Classes and Instances Nowadays, it is a common practice in optimization to choose a fixed set of problem instances in advance and to apply classical ANOVA or regression analysis. In many experimental studies a few problem instances πi(i= 1,2, . . . , q) are used and results of some runs of the algorithms αj(j= 1,2, . . . , h) on these instances are collected. The instances can be treated as blocks and all algorithms are run on each single instance. Results are grouped per instance πi. Analyses of these experiments shed some light on the performance of the algorithms on those specific instances. However, the interest of the researcher should not be just the performance of the algorithms on those specific instances chosen, but rather on the generalization of the results to the entire class Π. Generalizations 3 about the algorithm’s performance on new problem instances are difficult or impossible in this setting. Based on ideas from Chiarandini and Goegebeur [11], to overcome this difficulty, we propose the following approach: A small set of problem instances {πi∈Π|i= 1,2, . . . , q}is chosen at random from a large set, or class Π, of possible instances of the problem. Problem instances are considered as factor levels. However, this factor is of a different nature from the fixed algorithmic factors in the classical ANOVA setting. Indeed, the levels are chosen at random and the interest is not in these specific levels but in the problem class Π from which they are sampled. Therefore, the levels and the factor are random. Consequently, our results are not based on a limited, fixed number of problem instances. They are randomly drawn from an infinite set, which enables generalization. 2.2 Feature Extraction and Instance Generation A problem class Π can be generated in different manners. We will consider artificial and natural problem class generators. Artificially generated problems allow feature generation based on some predefined characteristics. They are basically theory driven, i.e., the researcher defines certain features such as linearity or multi modality. Based on these features, a model (formula) is constructed. By integrating parameters into this formula, many problem instances can be generated by parameter variation. We will exemplify this approach in the following paragraph. The second way, which will generate natural problem classes, uses a three-stage approach. First, the real-word system and its components are described. Then, features are extracted from a real-world system. Based on this feature set, a model is defined. Adding parameters to this model, new problem instances can be generated. There is also a third way to ”generate” test instances: if we are lucky, many data are available. In this case, we can sample a limited number of problem instances from the larger set of real-world data. The statistical analysis is similar for these three cases. 2.2.1 Artificial Test Functions Several problem instance generators have been proposed over the last years. For example, Gallagher and Yuan present a landscape test generator, which can be used to set up problem instances for continuous, bound-constrained optimization problems [13]. The Max-Set of Gaussian Landscape Generator (MSG) uses the maximum of mweighted Gaussian functions G(x) = max i∈1,2,...,m(wigi(x)), where g:Rn→Rdenotes an n-dimensional Gaussian function g(x) = exp −1 2(x−µ)Σ−1(x−µ)T (2π)n/2|Σ|1/2!1/n , 4 µis an n-dimensional vector of means, and Σ is an (n×n) covariance matrix. The mean of each Gaussian corresponds to an optimum on the landscape and the location of all optima is known. The global optimum is the one with the largest value. We will use the MSG problem instance generator in Sect. 5 to demonstrate our approach. 2.2.2 Natural Problem Classes This section exemplifies the three fundamental steps for generating real-world problem instances, namely 1. Describing the real-world system and its data 2. Feature extraction and model construction 3. Instance generation. We will illustrate this procedure by using the classic Box and Jenkins airline data [8]. These data contain the monthly totals of international airline passengers from 1949 to 1961. The feature extraction is based on methods from time-series analysis. Because of its simplicity the Holt-Winters method is popular in many application domains. It is able to adapt to changes in trends and seasonal patterns. The Holt-Winters prediction function requires the estimation of three parameters, i.e., α,βand γ, which can be estimated from original timeseries data. Their optimal values are determined by minimizing the squared one-step prediction error. To generate new problem instances, these parameters can be slightly modified. Based on these modified values, the model is re-fitted. Finally, we can extract the new time series. One typical result from this instance generation is shown in Fig. 1. Bartz-Beielstein [2] describes this procedure in detail. To illustrate the wide applicability of this approach, we will list further realwork problem domains, which are subject of our current research. Smart Metering. The development of accurate forecasting methods for electrical energy consumption profiles is an important task. We consider time series collected from a manufacturing process. Each time series contains quarter-hourly samples of the energy consumption of a bakery. A detailed data description can be found in [3]. Water Industry. Canary is a software developed by the United States Environmental Protection Agency (US EPA) and Sandia National Laboratories. Its purpose is to detect events in the context of water contamination. An event is in this context defined as a certain time period where a contaminant significantly deteriorates the water quality. Distinguishing events from (i) background changes, (ii) maintenance and modification due to operation, and (iii) outliers is an essential task, which was implemented in the Canary software. Therefore, deviations are compared to regular patterns and short term changes. The corresponding data contains multivariate time-series data. It is a selection from a larger dataset shipped with the open source event-detection software CANARY developed by US EPA and Sandia National Laboratories [19]. 5 Time AirPassengers 1950 1952 1954 1956 1958 1960 100 200 300 400 500 600 Figure 1: Holt-Winters problem instance generator. The solid line represents the real data, the dotted line predictions from the Holt-Winters model and the fine dotted line modified predictions, respectively. Finance. The data are real-world data from intraday foreign exchange (FX) trading. The FX market is a financial market for trading currencies to enable international trade and investment. It is the largest and most liquid financial market in the world. Currencies can be traded via a wide variety of different financial instruments, ranging from simple spot trades over to highly complex derivatives. We are using three foreign exchange (currency rate) time series collected from Bloomberg. Each time series contains hourly samples of the change in currency exchange rate [12]. One typical goal in forecasting is the minimization of the forecast errors or the differences between real (observed) values, say yi, and predicted values, say ˆyi. This goal can be considered as an optimization problem. As stated in Sect. 2.2, the statistical analysis is similar for artificial and natural problem classes. Our goal can be stated as follows: For a given problem class Π, which can be artificial or natural, we are trying to determine if an optimization algorithm αor several algorithm instances αishow similar behavior on randomly selected problem instances πi∈Π. This question will be formulated as a statistical hypothesis. Based on the related statistical framework, we can determine confidence intervals for the performance of the algorithm on unseen problem instances. 6 3 Algorithm Features 3.1 Factors and Levels Evolutionary algorithms (EA) belong to the large class of bio-inspired search heuristics. They combine specific components, which may be qualitative, like the recombination operator or quantitative, like the population size. Our interest is in understanding the contribution of these components. In statistical terms, these components are called factors. The interest is in the effects of the specific levels chosen for these factors. Hence, we say that the levels and consequently the factors are fixed. Although modern search techniques like sequential parameter optimization or Pareto genetic programming [18] allow multi-objective performance measures (solution quality versus variability or description length), we restrict ourselves to analyze the effect of these factors on a univariate measure of performance. We will use the quality of the solutions returned by the algorithm at termination as the performance measure. 3.2 Example: Evolution Strategy Evolution strategies (ES) are prominent representatives of evolutionary algorithms, which includes genetic algorithms and genetic programming as well [17]. They can be classified as generic population-based metaheuristic optimization algorithms for global optimization that in some sense mimics the natural evolution. Evolution strategies are applied to hard real-valued optimization problems. Mutation is performed by adding a normally distributed random value to each vector component. The standard deviation of these random values is modified by self-adaptation. Evolution strategies can use a population of several solutions. Each solution is considered as as individual and consists of object and strategy variables. Object variables represent the position in the search space, whereas strategy variables store the step sizes, i.e., the standard deviations for the mutation. We are analyzing the ES basic variant, which has been proposed in [6]. Mutation means neighborhood-based movement in search space that includes the exploration of the ”outer space” currently not covered by a population, whereas recombination rearranges existing information and so focuses on the ”inner space”. Selection is meant to introduce a bias towards better fitness values. A concrete ES may contain specific mutation, recombination, or selection operators, or call them only with a certain probability, but the control flow is usually left unchanged. Each of the consecutive cycles is termed a generation. The control flow is shown in Fig. 2. Concerning the representation, it should be noted that most empiric studies are based on canonical forms as binary strings or real-valued vectors, whereas many real-world applications require specialized, problem dependent ones. Table 1 summarizes important ES parameters. This chapter presents two case studies. The first case study is based on a fixed ES parameter setting, whereas the second case study modifies the recombination operator for object variables. We are convinced that the applicability of the 7 mating selection recombination initialization and evaluation mutationevaluation test for termination environmental selection crossover replacement Figure 2: The evolutionary cycle, basic working scheme of all ES and EA. Terms common for describing evolution strategies are used, alternative terms are added below in blue. methods presented in this chapter goes far beyond the simplified case studies. Our main contribution is a framework, which allows conclusions that are not limited to a small number of problem instances but to problem classes. 4 Objective Functions We will use the following optimization framework: An ES is applied as a minimizer on the test function f(x). Formally speaking, let Sdenote some set, e.g., S⊆Rn. We are seeking for values f∗and x∗, such that minx∈Sf(x) with f∗= minx∈Sf(x) and x∗= arg min f(x). This approach can be extended in many ways. For example, if Sdenotes times-series data, then an optimization algorithm can be applied to minimize the empirical mean squared prediction error. Test problem instances will be drawn from Gallagher’s and Yuan’s MSG test function generator. The following parameters can be used to specify the MSG generator. •The number of Gaussian components m. •The mean vector µof each component. •The covariance matrix Σ of each component. •The weight of each component wi. •A maximum threshold t∈[0; 1] can be specified for local optima and the fitness value of the global optimum G∗. Local optima are randomly generated within [0; t×G∗]. 8 Table 2: ANOVA table for a one-factor fixed and random effects models Source Sum Degrees Mean EMS EMS of Variation of Squares of freedom Square Fixed Random Treatment SStreat q−1 MStreat σ2+rPq i=1 τ2 i q−1σ2+rσ2 τ Error SSerr q(r−1) MSerr σ2σ2 Total SStotal qr −1 and MSerr =SSerr q(r−1) =Pq i=1 Pr j=1(Yij −¯ Yi.)2 q(r−1) . It can be shown that E(MStreat) = σ2+rσ2 τand E(MSerr) = σ2,(9) cf. [15]. Therefore, the estimators of the variance components are ˆσ2= MSerr,(10) ˆ σ2 τ=MStreat −MSerr r.(11) The corresponding ANOVA table is shown in Table 2. Based on ANOVA calculations, we obtain with (10) an estimator of the first variance component ˆσ2=−0.4848257, and from (11), we obtain the second component ˆσ2 τ= 11.32854. The model variance can be determined as ˆσ2+ ˆσ2 τ= 10.84372. The mean µ=−12.05554 from (8) can be extracted. Finally, the p value in the ANOVA table is calculated as 0.7979083. Note, that we have obtained a negative variance. Since negative variances are not feasible, we can proceed by setting their values to zero and proceed with this modified values. A more elegant way is presented in the following. Restricted maximum likelihood. In some cases, the standard ANOVA, which was used in our example, produces a negative estimate of a variance component. This can be seen in (11): If MSerr >MStreat, negative values occur. By definition, variance components are positive. Methods, which always yield positive variance components have been developed. Here, we will use restricted maximum likelihood estimators (REML). The ANOVA method of variance component estimation, which is a method of moments procedure, and REML estimation may lead to different results. Output from an R-based analysis with the function lme from the package lme4 reads as follows (fSeed denotes the problem instance) [16]: Linear mixed model fit by REML Formula: yLog ~ 1 + (1 | fSeed) Data: samp.df 15 ●● ● ●●● ●● ●● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●● ●● ● ● ● ●● ● ●● ●● ●● ●● ●● ● ●●● ●● ● ●●● ● ●● ● ●●● ●● ●● ●●● ● ● ●● ● ●● ● ● ●●● ● ● ●●● ● ●● −2 −1 0 1 2 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Q−Q plot for residuals Theoretical Quantiles Sample Quantiles ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 −1 0 1 2 −5 0 5 10 Q−Q plot for residuals Theoretical Quantiles Sample Quantiles Figure 5: Left: Q-Q plot of the residuals for raw data. Left: Q-Q plot for the log-transformed responses. AIC BIC logLik deviance REMLdev 475.6 483.1 -234.8 469.3 469.6 Random effects: Groups Name Variance Std.Dev. fSeed (Intercept) 0.000 0.0000 Residual 10.893 3.3004 Number of obs: 90, groups: fSeed, 9 Fixed effects: Estimate Std. Error t value (Intercept) -12.0555 0.3479 -34.65 Compared to the ANOVA setting, different values for ˆσ2, ˆσ2 τ, and µwere obtained. However, the REML based analysis also shows that the variability in the response observations can be attributed to the variability of the algorithm. SAMP-3 Validation of the Model Assumptions Before performing hypothesis testing based on the models introduced in SAMP-2, the validity of the model assumptions has to be investigated. If the model is adequate, the residuals should exhibit no structure. Residuals are plotted against fitted values to check the assumption of homoscedasticity and quantile-quantile (Q-Q) plots are used to check if residuals meet the normality assumption. Quantile-quantile plots of the residuals is shown in Fig. 5 for the raw and the log-transformed responses. These plots provide a good way to compare the the distribution of a sample with a distribution. Large deviations from the line indicate non-normality of the sample data. These Q-Q plots indicate that a log transformation of the response might be useful in our setting. SAMP-4 Hypothesis Testing Testing hypotheses about individual treatments (instances) is useless, because the problem instances πiare here consid16 ered as samples from some larger population of instances Π. We test hypotheses about the variance component σ2 τ, i.e., the null hypothesis H0:σ2 τ= 0 versus H1:σ2 τ>0.(12) Under H0, the algorithm performance is identical on all problem instances (”all treatments are identical”), i.e., rσ2 τis very small. Based on (9), we conclude that E(MStreat) = σ2+rσ2 τand E(MSerr) = σ2are similar. Under the alternative, variability exists between treatments. Standard analysis shows that SSerr/σ2is distributed as chi-square with q(r−1) degrees of freedom. Let Fu,v denote the Fdistribution with unumerator and vdenominator degrees of freedom. Under H0, the ratio F0= SStreat q−1 SSerr q(r−1) =MStreat MSerr is distributed as Fq−1,q(r−1). To test hypotheses in (8), we require that τ1, . . . , τq are i.i.d. N(0, σ2 τ), εij,i= 1, . . . , q,j= 1, . . . , r, are i.i.d. N(0, σ2), and all τi and εij are independent of each other. These considerations lead to the decision rule to reject H0at the significance level αif f0> F(1 −α;q−1, q(r−1)),(13) where f0is the realization of F0from the observed data. An intuitive motivation for the form of statistic F0can be obtained from the expected mean squares. Under H0both MStreat and MSerr estimate σ2in an unbiased way, and F0can be expected to be close to one. On the other hand, large values of F0give evidence against H0. Regarding the SAMP case, we obtain the following values: Based on (9) and (13), we can determine the Fstatistic and the pvalue. We get MStreat = MSerr = 10.89275 and f0= 1, which results a large pvalue: 0.4426363. The null hypothesis H0:σ2 τ= 0 from (12) can not be rejected, i.e., we conclude that there is no instance effect. A similar conclusion was obtained from the ANOVA method of variance component estimation as introduced in Table 2. SAMP-5 Confidence Intervals and Prediction An unbiased estimator of the overall mean µis ˆµ= ¯y·· = q X i=1 r X j=1 yij/(qr). Its variance is given by V(¯y··)=V  q X i=1 r X j=1 yij/(qr) =rσ2 τ+σ2 qr . With (9) and (10), we obtain an estimator of the variance of the overall mean µas ˆ V(¯y··) = MStreat/qr. 17 Since ¯ Y·· −µ pMStreat/qr ∼tq(r−1), the confidence limits for µcan be derived as ¯y·· ±t1−α/2;q(r−1)pMStreat/qr. (14) We conclude the SAMP case study with prediction of the algorithm’s performance on a new instance from the same class. Based on (14), we obtain the following 95% confidence interval: [2.6773e−06; 1.262e−05]. Again, confidence intervals from the REML and ANOVA methods are very similar. Summarizing, we can conclude that the ES performs similar on instances from ΠMSG, which were generated with Eq. 2. 5.3 MAMP: Multiple Algorithms, Multiple Problems In the MAMP case study, fixed effects are included in the conditional structure of (6), which leads to a mixed model. Instead of one fixed algorithm as in the SAMP case, we consider either several algorithms or algorithms with several parameters. Both situations can be treated while considering algorithms as levels of a fixed factor, whereas problem instances are drawn randomly from the population of instances ΠMSG. MAMP-1 Algorithm and Problem Instances MAMP-2 ANOVA and REML Model Building MAMP-3 Validation of the Model Assumptions MAMP-4 Hypothesis Testing a) Random effects b) Fixed effects MAMP-5 Confidence Intervals and Prediction MAMP-1 Algorithm and Problem Instances We aim at comparing the performance of the ES with different recombination operators over an instance class. More precisely, we have four ES instances using recombination operators {1,2,3,4}and nine instances randomly sampled from the class ΠMSG as illustrated in Fig. 3. Each run is repeated ten times. In this study 4 ×9×10 = 360 data were used. We are interested in the following questions: •Is there an instance effect? •Do the mean performances of the ES with different recombination operators differ? •Do the instance-algorithm interactions contribute to the variability of the response? A first visual inspection, which plots the performance of the algorithm within each problem instance, is shown in Fig. 6. In eight of the nine instances the linear regression line does have a negative slope and the intercepts do not differ very much. This indicates that there is no significant interaction between the fixed and the random factors. 18 objreco y −10 −5 1 2 3 4 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 124 ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 125 1234 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 126 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 127 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 128 −10 −5 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 129 −10 −5 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● 130 1234 ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 131 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 132 Figure 6: Four algorithms (ES with modified recombination operators) on nine test problem instances. Each panel represents one problem instance and problem instances are labeled from 124 to 130. Performance is plotted against the level of the recombination operator. 19 MAMP-2 ANOVA and REML Model Building The variability in the performance measure can be decomposed according to the following mixedeffects ANOVA model: Yijk =µ+αj+τi+γij +εijk,(15) where µis an overall performance level common to all observations, αjis a fixed effect due to the algorithm j,τiis a random effect associated with instance i, γij is a random interaction between instance iand algorithm j, and εijk is a random error for replication kof algorithm jon instance i. We assume that the αj’s are fixed effects such that Ph j=1 αj= 0 and that the random elements are τiare i.i.d. N(0, σ2 τ), γij are i.i.d. N(0, σ2 γ), εijk are i.i.d. N(0, σ2), and τi,γij and εijk are mutually independent random variables. Similar to (6) the conditional distribution of the performance measure given the instance and the instance–algorithm interaction is given by Yijk|τi, γij ∼ N(µ+αj+τi+γij, σ2),(16) with i= 1, . . . , q, j = 1, . . . , h, and k= 1, . . . , r. The marginal model reads (after integrating out the random effects τiand γij): Yijk ∼ N (µ+αj, σ2+σ2 τ+σ2 γ).(17) Based on these statistical assumptions, hypothesis tests can be performed about fixed and random factor effects. Using the mixed model (16), we are interested in testing whether there is a difference between the factor level means µ+αj (j= 1, . . . , h). The hypotheses for testing the fixed effects can be formulated as H0:αi= 0 ∀iagainst H1:∃αj6= 0 (18) Regarding random effects, tests about particular levels are useless. This is similar to the random-effects model (8). Again, we perform tests about the variance components σ2 τand σ2 γinstead. These can be formulated as follows: H0:σ2 τ= 0,and H0:σ2 γ= 0, H1:σ2 τ>0, H1:σ2 γ>0,(19) respectively. If all treatment (problem instances) combinations have the same number of observations, i.e., if the design is balanced, the test statistics for these hypotheses are ratios of mean squares that are chosen such that the expected mean squares of the numerator differs from the expected mean squares of the denominator only by the variance components of the random factor under test. Chiarandini and Goegebeur [11] present the resulting analysis of variance, which is shown in Table 3. ANOVA Model Building. The ANOVA table for the experiments from the MAMP case study is shown in Table 4. Equating the observed mean squares 20 Table 3: Expected mean squares and consequent appropriate test statistics for a mixed two-factor model with hfixed factors, qrandom factors, and rrepeats. From [11]. Mean Expected Test Effects squares df mean squares statistics Fixed factor MSA h−1σ2+rσ2 γ+rq Ph j=1 α2 j h−1MSA/MSAB Random factor MSB q−1σ2+rσ2 γ+rhσ2 τMSB/MSAB Interaction MSAB (h−1)(q−1) σ2+rσ2 γMSAB/MSE Error MSE hq(r−1) σ2 Table 4: ANOVA for the MAMP case Mean squares Factors Df Sum Sq Mean Sq F value Pr(>F) MSA objreco 3 154.59 51.53 11.05 0.0000 MSB fSeed 8 251.79 31.47 6.75 0.0000 MSAB objreco:fSeed 24 185.60 7.73 1.66 0.0288 MSE Residuals 324 1511.27 4.66 21 in the lines of the ANOVA table to their expected values and solving for the variance components leads to the following equations[15]: ˆσ2 τ=MSB −MSAB hr = 0.593502 ˆσ2 γ=MSAB −MSE r= 0.306907 ˆσ2=MSE = 4.664423 Next, we will compare these results to the REML based analysis of the mixed model. REML Model Building. We have specified sum contrasts instead of the default treatment contrasts used in lmer(). Again, fSeed represents the problem instance, whereas the algorithm instance αj,j= 1,...,4, is represented by objreco. Linear mixed model fit by REML Formula: yLog ~ objreco + (1 | fSeed) + (1 | fSeed:objreco) Random effects: Groups Name Variance Std.Dev. fSeed:objreco (Intercept) 0.30691 0.55399 fSeed (Intercept) 0.59351 0.77039 Residual 4.66442 2.15973 Number of obs: 360, groups: fSeed:objreco, 36; fSeed, 9 Fixed effects: Estimate Std. Error t value (Intercept) -6.0222 0.2956 -20.370 objreco1 0.6176 0.2539 2.433 objreco2 0.6918 0.2539 2.725 objreco3 -0.6671 0.2539 -2.628 As can be seen from the Random effects section of the REML model output, the estimated variances for the problem instance and the instance-interaction random effects are ˆσ2 τ= 0.59351 and ˆσ2 γ= 0.30691, respectively. The Random effects section presents the estimates of the fixed effects model parameters, i.e., objreco. MAMP-3 Validation of the Model Assumptions Again, the check of the diagnostic plots (Fig. 7) reveals that a log transformation of the response improves the model adequacy. 22 ●●● ●● ● ●●● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ●●● ● ● ● ● ●● ● ● ●● ● ●● ● ● ●● ●●● ● ●● ●● ● ●●● ● ● ● ● ● ●● ● ● ● ●●● ● ●● ●●● ● ● ● ●●● ●● ● ● ●● ●● ● ●● ● ●● ● ● ● ● ● ● ●●● ●● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ●●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ●● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ●● ● ● ●●● ●● ●● ● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ●●● ● ●●● ● ● ●● ●●● ● ●● ● ●● ●● ●● ● ● ●● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ●●● ●●● ●● ● ● −3 −2 −1 0 1 2 3 0.0 0.1 0.2 0.3 Q−Q plot for residuals Theoretical Quantiles Sample Quantiles ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −3 −2 −1 0 1 2 3 −6 −4 −2 0 2 4 Q−Q plot for residuals Theoretical Quantiles Sample Quantiles Figure 7: Left: Q-Q plot of the residuals for raw data. Left: Q-Q plot for the log-transformed responses. MAMP-4a Hypothesis Testing: Random Effects We will consider random effects first. Regarding problem instances, test about levels are meaningless. Hence, we perform tests about the variance components σ2 τand σ2 γ, which were presented in (19). First, we are testing the null hypothesis, which states that the components of the random effects are zero. Based on the ANOVA from Table 3, we obtain the values for the MAMP case that are shown in Table 4. The values reveal that there are main factor effects (fixed and random), but no significant interaction effects. Alternatively, we can compute the likelihood ratios of models with and without the factors under observation. Data: mamp.df Models: mamp.lmer2: yLog ~ objreco + (1 | fSeed) mamp.lmer3: yLog ~ objreco + (1 | fSeed) + (1 | fSeed:objreco) Df AIC BIC logLik mamp.lmer2 6 1616.7 1640.0 -802.35 mamp.lmer3 7 1616.6 1643.8 -801.31 Chisq Chi Df Pr(>Chisq) 2.0929 1 0.148 These tests indicate that there are also no significant instance-algorithm interactions. Additional likelihood-ratio test show that the fixed factor and random factor effects are significant. MAMP-4b Hypothesis Testing: Fixed Factor Effects Regarding fixed factors, we are interested in testing for differences in the factor level means µ+αi. These tests were formulated in (18), i.e., we are testing H0: all αi are equal to 0 versus H1: at least one αj6= 0. Here, we are using the test 23 statistic from [15, p. 523] for testing that the means of the fixed factor effects are equal. The appropriate test statistic for testing that the means of the fixed factor effects are equal, i.e., H0is true, is F0=MSA MSAB =154.59/3 185.6/24 = 6.663362, with values taken from Table 4. The reference distribution is Fn−1,(n−1)(q−1). We calculate the pvalue for the test on the fixed-effect term. The obtained p value is 0.002, hence the results collected indicate that the factor recombination (objreco) has a statistically significant impact on the performance of the algorithm. Using sum of contrasts implies that Pαj= 0. The point estimates for the mean algorithm performance with the jth fixed factor setting can be obtained by µ·j=µ+αj. The fixed factor effects can be estimated in the mixed model as ˆµ=y... ˆαj=yj. −y..., which results in the following estimates: ˆα1= 0.6175519, ˆα2= 0.6918047, ˆα3=−0.6671266, and ˆα4=−0.6423659. The same estimates were obtained with the REML analysis as can be seen from the REML model output on page 22. The corresponding fixed effects are shown in the Fixed effects section of the REML output. For example, we obtain the following value: objreco1 = ˆα1= 0.6176. MAMP-5 Confidence Intervals and Prediction We generate paired comparisons plots, which are based on confidence intervals. The wrapper function intervals() from Chiarandini and Goegebeur [11] was used for visualizing these confidence intervals as shown in Fig. 8. When intervals overlap we conclude that there is no significant difference. Here, we can conclude that the recombination operators (1) and (2) show a similar performance, whereas performances between (3) and (2) are different. Intermediate recombination of the object variables, i.e., (3) and (4), results in a significant improvement of the performance. 6 Summary and Outlook In order to answer question (Q-1), we propose an approach to generate natural problem classes, which are based on real-world data. If no such data are available, artificial problem generators such as MSG can be used. Since our approach uses a model, say M, to generate new problem instances, one conceptual problem arises: This approach is not applicable, if the final goal is the determination of a model for the data, because Mis per definition the best model in this case and the search for good models will result in M. But there is a simple solution to this problem. In this case, the feature extraction and 24