micompr: An R Package for Multivariate Independent Comparison of Observations

The R package micompr implements a procedure for assessing if two or more multivariate samples are drawn from the same distribution. The procedure uses principal component analysis to convert multivariate observations into a set of linearly uncorrelated statistical measures, which are then compared using a number of statistical methods. This technique is independent of the distributional properties of samples and automatically selects features that best explain their differences. The procedure is appropriate for comparing samples of time series, images, spectrometric measures or similar high-dimension multivariate observations.


Introduction
The aim of this paper is to present the micompr package for R (R Core Team, 2015), which implements a procedure for comparing multivariate samples associated with different factor levels or groups. The research goal is to differentiate among pre-specified, well-defined classes or groups of sampling entities generating highly multivariate observations in which the dimensions or dependent variables are correlated, and to test for significant differences among groups. The procedure uses principal component analysis (PCA) (Jolliffe, 2002) to convert multivariate observations into a set of linearly uncorrelated statistical measures, which are then compared using a number of statistical methods, such as hypothesis tests and score plots.
This technique has several desirable attributes: a) it automatically selects observation features that best explain sample differences; b) it does not depend on the distributional properties of samples; and, c) it simplifies the researchers' work, as it can be used directly on multivariate observations. The procedure is appropriate for comparing samples of multivariate observations with highly correlated and similar scale dimensions, such as time series, images or spectrometric measures. However, the micompr package goes one step further by also accommodating the simultaneous comparison of multiple observation types, i.e., multiple outputs from a given "system". In this context, a "system" can be defined as an abstract entity capable of generating one or more stochastic data streams, i.e., outputs. Thus, micompr can determine if two or more instances of such a "system" display the same behavior by comparing observations of their outputs.
The remainder of this paper is organized as follows. First, in Section 2, commonly used techniques for differentiating samples of multivariate observations are discussed. The methodology employed by micompr is described in Section 3. Section 4 introduces the software and its architecture, namely the available objects and functions. Several concrete application examples, and how the "system"-output terminology fits each one, are presented in Section 5. The paper closes with Section 6, in which the overall approach and the R package are summarized.
1. Apply a MANOVA test to the samples, where each observation has q-dimensions, corresponding to the first q PCs (dimensions) such that these explain a user-defined minimum percentage of variance.

Apply a univariate test to observations in individual PCs. Possible tests include the t-test and the
Mann-Whitney U test for comparing two samples, or ANOVA and Kruskal-Wallis test, which are the respective parametric and non-parametric versions for comparing more than two samples.
The MANOVA test yields a single p-value from the simultaneous comparison of observations along multiple PCs. An equally succinct answer can be obtained with a univariate test using the Bonferroni correction or a similar method for handling p-values from multiple comparisons. However, both approaches will not prioritize dimensions, even though the first PCs can be more important for characterizing an output, as they explain more variance. In the univariate case one can prioritize PCs according to the explained variance using the weighted Bonferroni procedure (Rosenthal and Rubin, 1983).
Conclusions concerning whether samples are statistically similar can be drawn by analyzing the pvalues produced by the employed statistical tests, which should be below the typical 1% or 5% when samples are significantly different. In such case, less PCs should be required to explain the same percentage of variance than when, in the same context, no significant differences are found. The scatter plot of the first two PC dimensions can also provide visual, although subjective feedback on sample similarity.
While the procedure is most appropriate for comparing multivariate observations with highly correlated and similar scale dimensions, assessing the similarity of "systems" with multiple outputs of different scales is also possible. The simplest approach would be to apply the proposed method to samples of individual outputs, and analyze the results in a multiple comparison context. An alternative approach consists in concatenating, observation-wise, all outputs, centered and scaled to the same order of magnitude, thus reducing a "system" with k outputs to a "system" with one output. The proposed method would then be applied to samples composed of concatenated observations encompassing the existing outputs. This technique is described in detail by Fachada et al. (2017) in the context of comparing simulation outputs.

Overview
The micompr package for the R statistical computing environment implements the methodology proposed in Section 3. Here we describe version 1.0.1 of the package, which is available at https://cran. r-project.org/package=micompr. The development version is hosted at https://github.com/fakenmc/ micompr. The package is covered by the MIT license.
The micompr package is built upon two functions, cmpoutput and micomp. The former compares two or more samples of multivariate observations collected from one output. The latter is used for comparing multiple outputs and/or comparing outputs in multiple contexts. grpoutputs is a helper function for loading data from two or more set of files and preparing the data to be processed by the cmpoutput and/or micomp functions. assumptions is a generic function for assessing the assumptions of the parametric tests used in sample comparisons.

Architecture
micompr is structured according to the S3 object-oriented system. The cmpoutput, micomp and grpoutputs functions produce S3 objects with the same name. The package also provides the generic function assumptions, and two concrete implementations of methods for "cmpoutput" and "micomp" objects, which return objects of class "assumptions_cmpoutput" and "assumptions_micomp", respectively. All classes have method implementations of the common S3 generic functions print, summary and plot. Additionally, method implementations of the toLatex function, for producing user-configurable L A T E X tables with information about the performed comparisons, are provided for "cmpoutput" and "micomp" objects.

grpoutputs
This function groups outputs from sets of files containing multiple observations into samples. It returns a list of output matrices, ready to be processed by micomp. Alternatively, individual output matrices can be handled by cmpoutput. Separate files contain one multivariate observation of one or more outputs, one column per output, one row per dimension or variable. Each specified set of files is associated with a different sample. The function is also able to create an additional concatenated output, composed from the centered and scaled original outputs.
The plot method for "grpoutputs" objects shows k plots, one per output. Output observations are plotted on top of each other, with different samples colored distinctively. The summary method for "grpoutputs" objects returns a list containing two elements: a) the n × m dimensions of each output matrix; and, b) the sizes of individual samples. The print method for "grpoutputs" objects simply outputs the summary in a more adequate presentation format.

cmpoutput
The cmpoutput function is at the core of micompr. It compares two or more samples of multivariate observations using the technique described in Section 3. It accepts an output matrix, X (n×m) , with n observations and m variables or dimensions, a factor vector of length n, specifying the group associated with each observation, and a vector of explained variances with which to determine the number of PCs to use in the MANOVA test (alternatively, the number of PCs can also be directly specified). The function returns matrix T (n×r) of PCA scores and the p-values for the performed statistical tests, namely: a) a MANOVA test for each explained variance (or number of PCs); and, b) parametric (t-test or ANOVA) and non-parametric (Mann-Whitney or Kruskal-Wallis) univariate tests for each PC. Regarding the latter, the function also returns p-values adjusted with the weighted Bonferroni correction, using the percentages of explained variance by PC as weights.
The plot implementation for "cmpoutput" objects shows six sub-plots, namely a scatter plot with the PC1 vs. PC2 scores and five bar plots. The horizontal scale of the latter consists of the r PCs, and the vertical bars represent the explained variance (one plot) or univariate parametric and non-parametric pvalues, before and after weighted Bonferroni correction (four plots). The summary method for "cmpoutput" objects returns a list with the following items: a) percentage of variance explained by each PC; b) pvalues of the MANOVA test or tests; c) p-values of the parametric test, per PC, before and after weighted Bonferroni correction; d) p-values of the non-parametric test, per PC, before and after weighted Bonferroni correction; and, e) name of the parametric and non-parametric univariate tests employed (either t-test and Mann-Whitney U test for comparing two samples, or ANOVA and Kruskal-Wallis for more than two samples). The print method for "cmpoutput" objects shows the information provided by the summary implementation, but the p-values of the univariate tests are only shown for the first PC.

micomp
The micomp function performs one or more comparisons of multiple outputs, invoking cmpoutput for each comparison/output combination. It accepts a list of comparisons, where individual comparisons can have one of two configurations: a) a vector of folders and a vector of file sets containing data in the format required by grpoutputs, where each file set corresponds to a different sample; and, b) a "grpoutputs" object, passed directly. The returned objects, of class "micomp", are basically two-dimensional lists of "cmpoutput" instances, with rows associated with individual outputs, and columns with separate comparisons.
The plot method for "micomp" objects shows the PC1 vs. PC2 score plots for each comparison/output combination. The summary implementation for "micomp" objects returns a list of comparisons, each one containing a a × k matrix of p-values or number of PCs, associated with a ≥ 6 measures and k outputs. Four rows represent the p-values of the parametric and non-parametric univariate tests for the first PC, before and after weighted Bonferroni correction. The remaining pairs of rows are associated with the MANOVA test for a given percentage of variance to explain. One row shows the p-values, and the other displays the number of PCs required to explain the specified percentage of variance for the given output. As with other micompr objects, the print method for "micomp" objects also shows the summary with a better presentation.

assumptions
assumptions is a generic function which performs a number of statistical tests concerning the assumptions of the parametric tests performed by the package functions. Implementations of this generic function exist for "cmpoutput" and "micomp" objects. The former method returns objects of class "assumptions_cmpoutput" containing results of the assumptions tests for a single output comparison. The latter returns a two-dimensional list of "assumptions_cmpoutput" objects, with rows associated to individual outputs, and columns to separate comparisons. These objects are tagged with the "assumptions_micomp" class attribute.
The following assumptions are checked: a) observations are normally distributed within each sample along individual PCs (Shapiro-Wilk test); b) observations follow a multivariate normal distribution within each sample for all PCs used in MANOVA (Royston test); c) samples have homogeneous variance along individual PCs (Bartlett test); and, d) samples have homogeneous covariance matrices for all PCs used in MANOVA (Box's M test). Assumptions a) and c) should be verified for the parametric test applied to each PC, while assumptions b) and d) should be verified for individual MANOVA tests performed for each variance to explain (or, alternatively, for each specified number of PCs).
The plot implementations for classes "assumptions_cmpoutput" and "assumptions_micomp" display a number of bar plots for the p-values of the performed tests. These are more detailed for "assumptions_ cmpoutput" objects, showing the p-values of the univariate test for all PCs. For "assumptions_micomp" objects, one bar plot is shown per output/comparison combination, but in the case of the univariate tests only the p-values of the first PC are shown. Implementations of summary return a list of tabular data containing the p-values of the assumption tests. The summary method for "assumptions_cmpoutput" objects returns a list with two matrices of p-values, one for the MANOVA tests, another for the univariate tests. The summary method for "assumptions_micomp" objects follows the approach taken by the summary method for "micomp" objects, returning a list of p-value matrices, one matrix per comparison. Rows of individual matrices correspond to the assumptions tests, and columns to outputs. The print methods for "assumptions_cmpoutput" objects and for "assumptions_micomp" objects again show the summary information in a printable format.

toLatex methods for "cmpoutput" and "micomp" objects
These methods are implementations of the toLatex generic function, and convert "cmpoutput" and "micomp" objects to character vectors representing L A T E X tables. The generated tables are configurable via function arguments, with sensible defaults. Tables can present the following data for each output/comparison combination: a) number of principal components required to explain a user-specified percentage of variance; b) MANOVA p-value for a user-specified percentage of variance to explain or number of PCs; c) parametric test p-value for a given PC, before and/or after weighted Bonferroni correction; d) nonparametric test p-value for a given PC, before and/or after weighted Bonferroni correction; e) variance explained by a specific PC; and, f) a score plot with the output projection on the first two PCs.

Other functions
The micompr package is bundled with additional functions whose purpose is to aid the main package methods do their job. However, some of these may be useful in other contexts.
The concat_outputs function concatenates outputs collected from multiple observations. It accepts two arguments, namely a list of output matrices, and the centering and scaling method. Several centering and scaling methods, such as "range", "iqrange", "vast" or "pareto" (Berg et al., 2006), are recognized in the second argument. The function returns an n × p matrix of n observations with length p, which is the sum of individual output lengths. Lower-level centering and scaling of individual outputs is performed by the centerscale function, which accepts a numeric vector and returns a new vector, centered and scaled with the specified method.
The pvalf generic function formats p-values for L A T E X. A concrete default implementation is used by the micompr toLatex implementations. This implementation underlines and double-underlines p-values lower than 0.05 and 0.01, respectively, although these limits are configurable, and underlining can be turned off by setting both limits to zero. It is also possible to specify a limit below which p-values are capped. For example, if this limit is set to 1 × 10 −5 , a p-value equal to 1 × 10 −6 would be displayed as "< 1e −5 ". The default method of pvalf will format p-values lower than 5 × 10 −4 using scientific E notation, which is more compact and thus a better fit for tables. p-values between 5 × 10 −4 and 1 are formatted using regular decimal notation with three decimal places. This aspect is not configurable. However, another implementation of pvalf can be passed to the micompr toLatex implementations if different formatting is desired.
Simple Tik Z 2D scatter plots, as the ones produced by the micompr toLatex implementations, can be generated with the tikzscat function. The function accepts the data to plot, an n x 2 numeric matrix, of n observations and 2 dimensions, and a factor vector specifying the levels or groups associated with each observation. Several plot characteristics, such as mark types, scale and axes color, are configurable via function arguments. tikzscat returns a string containing the Tik Z figure code for plotting the specified data.

Included data
The package includes test data produced by several implementations of the Predator-Prey for High Performance Computing (PPHPC) simulation model (Fachada et al., 2015). The data is provided in rdata format, and is readily available on loading the package. The same data is also provided in TSV format. This is a limited subset of the complete data, and is included for package testing and exemplification purposes. The example discussed in Section 5.1 uses a superset of this data, which is available for public download, but could not be included with the package due to its large size.

Dependencies
micompr has a number of optional dependencies, not required for package installation and for using most of its functionality. The biotools (da Silva, 2016) and MVN (Korkmaz et al., 2014) packages are required by the assumptions functions, providing the statistical tests for assessing MANOVA and t-test assumptions. If these functions are invoked without the presence of the specified packages, they will inform the user of that fact, and terminate cleanly. The testthat (Wickham, 2011), knitr (Xie, 2015) and roxygen2 (Wickham et al., 2015) packages are required for package development. The deseasonalize package (McLeod and Gweon, 2013) is required for building one of the vignettes.

Examples
In this section we discuss four concrete application examples for the micompr package. The complete scripts used in these examples are available at https://github.com/fakenmc/micompr-examples.

Simulation model with multiple outputs
The replication of a simulation model in a new context highlights differences between the conceptual and implemented models, as well as inconsistencies in the conceptual model specification (Edmonds and Hales, 2003), promoting model verification, model validation (Wilensky and Rand, 2007), and model credibility (Thiele and Grimm, 2015). Some argue that a simulation model is untrustworthy until it has been successfully replicated (Edmonds and Hales, 2003;David, 2013). Model parallelization is a an illustrative example of the importance of replication, as it is often required for simulating large models in practical time frames (Fachada et al., 2016). By definition, model parallelization implies a number of changes, or even full reimplementation, of the original model, such that a robust comparison methodology, as provided by the micompr package, is required in order to make sure a parallelized model faithfully reproduces the behavior of the original serial model.
PPHPC is a reference model for studying and evaluating implementation strategies for spatial agentbased models, capturing important characteristics such as agent movement and local agent interactions (Fachada et al., 2015). The model describes a prototypical predator-prey system, and has six outputs, namely prey population, P s , predator population, P w , cell-bound food quantity, P c , mean prey energy, E s , mean predator energy, E w , and mean cell-bound food levels, C. Since outputs are collected once per iteration, each simulation run yields six time series, associated with the individual outputs. With several open source implementations publicly available, the model provides a good test case for multivariate comparison purposes.
Here we show the main comparison cases discussed in a previous article (Fachada et al., 2017), in which the model implementations are parameterized with size 400 and parameter set 1 (Fachada et al., 2015). A canonical PPHPC realization, implemented in NetLogo (Wilensky, 1999), is compared with three configurations of a parallel Java implementation (Fachada et al., 2016). The NetLogo implementation and the first Java configuration follow the PPHPC conceptual model and the specified parameters. The second Java configuration disables agent shuffling prior to agent actions, which is explicitly mandated in the conceptual model description. The third Java configuration performs a minimal change in one of the parameters specified by parameter set 1. As such, we define three comparison cases: Case I Compare the NetLogo implementation with the first Java configuration. These should yield distributionally equivalent results.
Case II Compare the NetLogo implementation with the second Java configuration. A small misalignment is to be expected.
Case III Compare the NetLogo implementation with the third Java configuration. There should be a mismatch in the outputs.
Independent samples of the six model outputs were obtained from n = 30 replications for each implementation or configuration, in a total of 4n = 120 runs. Each replication r = 1, . . . , 4n was performed with a PRNG seed obtained by taking the MD5 checksum of r, guaranteeing independence between seeds, and consequently, between replications. The data generated by this computational experiment, as well as the scripts used to set up the experiment, are made available to other researchers at https://zenodo.org/record/46848.
The following script performs these comparisons. Note that the concat = TRUE option of the micomp function specifies that an additional concatenated output, A, should be generated from the original outputs and analyzed in a similar fashion. The dir_data variable specifies the location of the dataset.

Monthly sunspots
This example uses the monthly sunspot data (WDC, 2016), included with R, which contains the monthly numbers of sunspots from 1749 to the present day. The solar cycle is an approximate 11-year period of changes in the number of sunspots and other associated phenomena. Thus, we divide the data into 11-year (132-month) periods, and consider each period to be an observation. In practice this is an oversimplification, since the cycles can be a bit longer or shorter than 11 years.
Given the data, we define two samples of 10 observations each, over a period of 110 years or 1320 months. The first sample includes solar cycles from 1749 to 1859, while the second encompasses cycles from 1902 to 2012. We can now ask the following question: were the solar cycles during the 1749-1859 interval significantly different from the more recent observations? The following code compares observations from the two time intervals, and attempts to provide an answer: # Factor vector, two levels: # a) ten 11-year cycles from 1749 to 1859 # b) ten 11-year cycles from 1902 to 2012 groups <-factor(c(rep("A", 10), rep("B", 10))) # Compare the two groups, use 9 PCs for MANOVA cmp <-cmpoutput("SunSpots", 9, m, groups) The cmp object can now be analyzed:

> cmp
Output name: SunSpots Number of PCs which explain 85.0 P-Value for MANOVA along 9 dimensions: 3.40755e-06 P-Value for t-test (1st PC): 2.713985e-06 P-Value for Mann-Whitney U test (1st PC): 4.330035e-05 Adjusted p-Value for t-test (1st PC): 6.513579e-06 Adjusted p-Value for Mann-Whitney U test (1st PC): 0.0001039211 The MANOVA p-value is significant, as well as the t-test and Mann-Whitney PC1 p-values, before and after weighted Bonferroni correction. As such, it is possible to conclude that solar cycles from 1749 to 1859 were significantly different from cycles between 1902 and 2012. However, is the data in accordance with the assumptions for the MANOVA and t-test? This can be checked with the assumptions function: Only Box's M test, which checks for homogeneity of variance-covariance matrices, is significant. However, this test is prone to false positives, and this assumption is not critical when samples are of the same size (Tabachnick and Fidell, 2013). Given this information, it seems plausible to consider the results provided by the parametric tests in our final decision, i.e., that there is in fact a significant difference between samples. A good way to visualize the overall results is to plot the "cmpoutput" object:

> plot(cmp)
This command generates the plots shown in Figure 1. The score plot shows the samples to be distinctly separated, and the variance explained by PC decreases abruptly from the first PC to the second. Univariate p-values for PC1 are visibly significant, though not very much for the remaining PCs.

Saugeen river flow
This example uses the Saugeen River daily flow data (Hipel and McLeod, 1994), included with the deseasonalize package. This data consists of a time series of the rivers' daily flow (m 3 /s) from 1915 to 1979. Considering one year as an observation, there are a total of 65 observations. We can, for example, define two samples of 30 observations each, with the first and last 30 years of records, and ask the following question: is there any statistical difference between the flow dynamics during the 1915-1944 and 1950-1979  # Factor vector, two levels: first 30 years and last 30 years groups <-factor(c(rep("A", 30), rep("B", 30))) # Compare cmp <-cmpoutput("SaugeenFlow", 0.9, saugdata, groups) The cmp object can now be analyzed: There are no significant p-values. As such, it is not possible to conclude that the river flow dynamic during the first 30 years of measurements is statistically different from the last 30 years. The MANOVA assumptions are not verified, as shown by be assumptions function: The t-test assumptions mostly hold, in spite of the PC1 normality test for sample A (first 30 years) being significant at the α = 0.05 level. Nonetheless, the U test p-value provides a similar conclusion. Plotting the cmp object provides another perspective: The previous command will produce plots in Figure 2. The PC1 vs. PC2 score plot does not show any clear sample separation and the decrease in explained variance from the first PC to the second is considerably less abrupt than what was observed for the sunspots example. Additionally, no t and U test p-values are significant for the 60 PCs after weighted Bonferroni correction, further reinforcing the conclusion that the yearly Saugeen river flow dynamic was similar during the compared periods of time.

PH 2 database of dermoscopic images
In this example we use the tools provided by the micompr package to study the PH 2 database of dermoscopic images (Mendonça et al., 2013). This image database contains a total of 200 dermoscopic images of melanocytic lesions, including, from benign to more serious, 80 common nevi, 80 atypical nevi, and 40 melanomas. These are 8-bit RGB color images, with a resolution of purportedly 768 × 560 pixels. We have found, however, that resolutions vary between 761 × 570 and 769 × 577. As such, we resized all images to 760 × 570 prior to our analysis. The goal is to verify if images of the three types of lesions form statistically distinguishable samples, i.e., this is not a classification exercise such as performed by Barata et al. (2014). Each image is considered an observation of three outputs, red, green and blue, corresponding to the respective color channels. The concatenation of all outputs, i.e., channels, provides a 4 th output. The three lesion samples are compared pairwise, as follows: 1v2 Common nevi and atypical nevi.
The following code reads the image dataset from disk and compares images grouped by lesion type. The imgfolder variable specifies the path containing the images (resized to 760×570), while the grpsfile variable specifies the path to the file containing the sample to which each image belongs.

Summary
In this paper we presented the R package micompr, which implements a procedure for comparing multivariate samples associated with different factor levels or groups. The package architecture and its core components were discussed and four examples were examined.