Analysis of Multivariate Data and Repeated Measures Designs with the R Package MANOVA.RM

The numerical availability of statistical inference methods for a modern and robust analysis of longitudinal- and multivariate data in factorial experiments is an essential element in research and education. While existing approaches that rely on specific distributional assumptions of the data (multivariate normality and/or characteristic covariance matrices) are implemented in statistical software packages, there is a need for user-friendly software that can be used for the analysis of data that do not fulfill the aforementioned assumptions and provide accurate p-value and confidence interval estimates. Therefore, newly developed statistical methods for the analysis of repeated measures designs and multivariate data that neither assume multivariate normality nor specific covariance matrices have been implemented in the freely available R-package MANOVA.RM. The package is equipped with a graphical user interface for plausible applications in academia and other educational purpose. Several motivating examples illustrate the application of the methods.


Introduction
Nowadays, a large amount of measurements are taken per experimental unit or subject in many experimental studies-requiring inferential methods from multivariate analysis in a unified way. Here we distinguish between two cases: 1. If the same quantity is measured under different treatment conditions or at different time points, a repeated measures (RM) design is present. Therein, observations are measured on the same scale and are combinable. This is also the case if the measuring instrument produces multiple responses, e.g., microarrays in bioinformatics.
2. If different quantities are measured on the same unit or subject, a multivariate analysis of variance (MANOVA) design is apparent. In such a situation, data is measured on different scales and not combinable (e.g., height and weight).
These two different definitions do not only lead to different questions of interest but also require different inference procedures as outlined below. In particular, the main difference between the two approaches is that in repeated measures designs comparisons between the response variables are meaningful. This means that also hypotheses regarding sub-plot or within subject factors (e.g., time) are of interest. On the other hand, MANOVA settings are usually designed to detect effects of the observed factors (and interactions thereof) on the multivariate outcome vectors, thus allowing -in contrast to multiple univariate ANOVA analyses-to evaluate the combined changes of the outcome variables with respect to the factor levels.
Despite their differences, MANOVA-and RM-type techniques share the same advantages over classical univariate endpoint-wise-ANOVA-type-analyses: • They provide joint inference and take the dependency across the endpoints into account, thus leading to possibly larger power to detect underlying effects.
• They allow for testing of additional factorial structures and • can easily be equipped with a closed testing procedure for subsequently detecting local effects in specific components, i.e. to perform post-hoc analyses.
Focusing on metric data and mean-based procedures, MANOVA and RM models are typically inferred by means of "classical" procedures such as Wilks' Lambda, Lawley-Hotelling, Roy's largest root (Davis, 2002;Johnson and Wichern, 2007;Anderson, 2001) or (generalized) linear mixed models with generalized estimating equations. For the classical one-way layout, these methods are implemented in R within the manova function in the stats package, where one can choose between the options 'Pillai', 'Wilks', 'Hotelling-Lawley' and 'Roy'. Nonparametric rank-based methods for null hypotheses formulated in distribution functions are implemented within the packages npmv for oneand two-way MANOVA (Burchett et al., 2017) and nparLD for several repeated measures designs (Noguchi et al., 2012). In case of fixed block effects, the GFD package (Friedrich et al., 2017b), which implements a permutation Wald-type test in the univariate setting, can also be used.
(Generalized) linear mixed models are implemented in the lm and the glm function (package stats) for univariate data as well as in the SCGLR package for Generalized Linear Model estimation in the context of multivariate data (Cornu et al., 2018). The Anova and Manova function in the car package (Fox and Weisberg, 2011) calculate type-II and type-III analysis-of-variance tables for objects produced by, e.g., lm, glm or manova in the univariate and multivariate context, respectively. In the MANOVA context, repeated measures designs can be included as well.
Most of these procedures, however, rely on specific distributional assumptions (such as multivariate normality) and/or specific covariance or correlation structures (e.g., homogeneity between groups or, for RM, compound symmetry; possibly implying equal correlation between measurements) which may often not be justifiable in real data. In particular, with decreasing sample sizes and increasing dimensions, such presumptions are almost impossible to verify in practice and may lead to inflated type-I-errors, cf. Vallejo et al. (2001); Lix and Keselman (2004) ;Vallejo Seco et al. (2007); Livacic-Rojas et al. (2010). To this end, several alternative procedures have been developed that tackle the above problems and have been compared in extensive simulation studies, see amongst others Brunner (2001); Lix and Lloyd (2007); Gupta et al. (2008) (2018) and the references cited therein. Here, we focus on nonparametric statistical methods that are valid in the multivariate Behrens-Fisher situation-equal covariance matrices across the groups are not assumed-and provide accurate inferential results in terms of p value estimates and confidence intervals for the parameters of interest. In particular, we implemented bootstrap-and permutation-based approaches to approximate the distribution of the test statistics in a robust way. Simulation studies comparing these approaches to the traditional methods mentioned above can, e. g. be found in the main papers and the supplements of Friedrich et al. (2017a) and Bathke et al. (2018).
More precisely, we focus on nonparametric methods for testing main and interaction effects of fixed factors in repeated measures designs and multivariate data. In particular, general Wald-type test statistics (for MANOVA and RM), ANOVA-type statistics (for RM) and modified ANOVA-type tests (for MANOVA) are implemented in MANOVA.RM (Friedrich et al., 2019) because • they can be used to test hypotheses in various factorial designs in a flexible way, • their sampling distribution can be approximated by resampling techniques, even allowing their application for small sample sizes, • and they are appropriate methods in the Behrens-Fisher situation.
To make the methods freely accessible we have provided the R package MANOVA.RM for routine statistical analyses. It is available from the R Archive at https://CRAN.R-project.org/package=MANOVA.RM The main functions RM (for RM designs) and MANOVA (for MANOVA designs) are developed in style of the well known ANOVA functions lm or aov. Its user-friendly application not only provides the p values and test statistics of interest but also a descriptive overview together with component-wise two-sided confidence intervals. Moreover, the MANOVA function even allows for an easy calculation and confidence ellipsoid plots for specified multivariate contrasts as described in .
Specifically, for testing multivariate main-and interaction effects in one-, two-and higher-way MANOVA models, the MANOVA function provides • the Wald-type statistic (WTS) proposed by Konietschke et al. (2015) using a parametric bootstrap, a wild bootstrap or its asymptotic χ 2 -distribution for p value computations, and • the modified ANOVA-type statistic (MATS) proposed by Friedrich and Pauly (2018) using a parametric or wild bootstrap procedure for p value computations.
In addition to multivariate group-wise effects, the RM function also allows to test hypotheses formulated across within subject factors. The implemented test statistics are • the ANOVA-type statistic (ATS) using an F-approximation as considered in Brunner (2001) as well as a parametric and a wild bootstrap approach and • the Wald-type statistic (WTS) using the asymptotic χ 2 -distribution (Brunner, 2001), the permutation technique proposed in Friedrich et al. (2017a) as well as a parametric (Bathke et al., 2018) and a wild bootstrap approach for p value estimation.
The paper is organized as follows: In Section 2.2 the multivariate statistical model as well as the implemented inference procedures are described. The application of the R package MANOVA.RM is exemplified on several Repeated Measures and MANOVA Examples in Section 2.3. Finally, the paper closes with a discussion in Section 2.4.
Throughout the paper we use the subsequent notation from multivariate linear models: For a ∈ N we denote by P a = I a − 1 a J a the a-dimensional centering matrix, by I a the a-dimensional identity matrix and by J a the a × a matrix of 1's, i.e., J a = 1 a 1 a , where 1 a = (1, . . . , 1) is the a-dimensional column vector of 1's.

Statistical model and inference methods
For both the RM and the MANOVA design equipped with an arbitrary number of fixed factors, we consider the general linear model given by d-variate random vectors (1) Here, k = 1, . . . , n i denotes the experimental unit or subject in group i = 1, . . . , a. Note, that a higher-way factorial structure on the groups/between subject or within subject factors can be achieved by sub-indexing the indices i (group/between subject factors) or j (within subject factors) into i 1 , . . . , i p or j 1 , . . . , j q . In this model µ i = (µ i1 , . . . , µ id ) ∈ R d is the mean vector in group i = 1, . . . , a and for each fixed i it is assumed that the error terms ik , k = 1, . . . , n i , are independent and identically distributed d-variate random vectors with mean E( i1 ) = 0 and existing variances 0 < σ 2 ij = var(X ijk ) < ∞, j = 1, . . . , d. For the WTS-type procedures we additionally assume positive definite covariance matrices cov( i1 ) = V i > 0 and existing finite fourth moments E(|| i1 || 4 ) < ∞.
Within this framework, hypotheses for RM or MANOVA can be formulated by means of an adequate contrast hypothesis matrix H by Let X • = (X 1· , . . . , X a· ) denote the vector of pooled group means X i· = 1 n i ∑ n i k=1 X ik , i = 1, . . . , a and Σ N = N · diag{ V 1 /n 1 , . . . , V a /n a } the estimated covariance of √ NX • . Here, N = ∑ i n i and In this set-up Konietschke et al. (2015) propose a statistic of Wald-type (WTS) for testing H 0 , where T = H (H H ) + H, X = {X 11 , . . . , X an a }, and A + denotes the Moore-Penrose inverse of the matrix A. Since its asymptotic χ 2 rank(T) null distribution provides a poor finite sample approximation, they propose the following asymptotic model-based bootstrap approach: Given the data X let X ik ∼ N(0, V i ), i = 1, . . . , a, k = 1, . . . , n i , be independent random vectors that are used for recalculating the test statistic as T N = T N (X ), where X = {X 11 , . . . , X an a }. Denoting by c the corresponding (1 − α)-quantile of the (conditional) distribution of T N the test rejects H 0 if T N > c . The validity of this procedure (also named parametric bootstrap WTS) is proven in Konietschke et al. (2015).
This procedure is not only applicable for MANOVA but also for RM designs (Bathke et al., 2018). However, Friedrich et al. (2017a) proposed a more favourable technique for Repeated Measurements. It is based on an at first blush chaotic resampling method: Wild permutation of all pooled components without taking group membership or possible dependencies into account. Denoting the resulting permuted data set as X π their permutation test for RM models rejects H 0 if T N > c π . Here c π denotes the (1 − α)-quantile of the (conditional) distribution of the permutation version of the test statistic T π N = T N (X π ). As shown in extensive simulations in Friedrich et al. (2017a) and the corresponding supporting information this 'wild' permutation WTS method controls the type-I error rate very well. Note that this procedure is only applicable for RM due to the commensurate nature of their components. In MANOVA set-ups the permutation would stir different scalings making comparisons meaningless.
In addition to these WTS procedures two other statistics are considered as well. For RM the The R Journal Vol. 11/2, December 2019 ISSN 2073-4859 well-established ANOVA-type statistic (ATS) by Brunner (2001) is implemented together with the enhanced F-approximation of the statistic proposed in Brunner et al. (1997Brunner et al. ( , 2002 and implemented in the SAS PROC Mixed procedure. Although known to be rather conservative it has the advantage (over the WTS) of being applicable in case of eventually singular covariance matrices V i or V i since it waives the Moore-Penrose inverse involved in Equation 2.
Similar to the permuted WTS the ATS given in Equation 3 is only applicable for RM since it is not invariant under scale transformations (e.g., change of units) of the univariate components. To nevertheless provide a robust method for MANOVA settings which is also applicable in case of singular V i or V i ,  have recently proposed the novel MATS (modified ATS) Here, the involved diagonal matrix D N = ⊕ 1≤i≤a,1≤s≤d N σ 2 is /n i of the empirical variances σ 2 is of component s in group i, deduces an invariance under component-wise scale transformations of the MATS for null hypotheses as described in Section 2.2.1, i.e., of the form T = M ⊗ I d , see  for details. To obtain an accurate finite sample performance, it is also equipped with an asymptotic model based bootstrap approach. That is, MATS rejects H 0 if M N >c , wherec is the (1 − α)-quantile of the (conditional) distribution of the bootstrapped statistic M N = M N (X ). In addition, we implemented a wild bootstrap approach, which is based on multiplying the centered data vectors In the package, we implemented Rademacher distributed weights, i.e., random signs. Extensive simulations in Friedrich and Pauly (2018) not only confirm its applicability in case of singular covariance matrices but also disclose a very robust behaviour that even seems to be advantageous over the parametrically bootstrapped WTS of Konietschke et al. (2015). However, both procedures, as well as the 'usual' asymptotic WTS are displayed within the MANOVA functions. All of the aforementioned procedures are applicable in various factorial designs in a unified way, i.e., when more than one factor may impact the response. The specific models and the hypotheses being tested will be discussed in the next section.

Special designs and hypotheses
In order to provide a general overview of different statistical designs and layouts that can be analyzed with MANOVA.RM we exemplify few designs that occur frequently in practical applications and discuss the model, hypotheses and limitations. All of the methods implemented in MANOVA.RM are even applicable in higher-way layouts than being presented here; and the list should not be seen as the limited application of the package. The models are derived by sub-indexing the index i in Equation 1 in the following ways: and obtain the null hypothesis of 'no group' or 'factor A' effect as In case of a = 2 this includes the famous multivariate Behrens-Fisher problem as, e.g., analyzed in Yao (1965); Nel and Van der Merwe (1986); Christensen and Rencher (1997); Krishnamoorthy and Yu (2004) or Yanagihara and Yuan (2005).
• Crossed Two-Way (A × B): Splitting the index into two and writing The corresponding null hypotheses of no main effects in A or B and no interaction effect between A and B can be written as: Note that the interpretation of main effects is complicated by the presence of significant interaction effects and further analyses are necessary to determine the direction of the effects.
• Hierarchically nested Two-Way (B(A)): A fixed subcategory B within factor A can be introduced via the model Here, the hypotheses of no main effect A or no sub-category main effect B can be written as .., µ 1b 1 , µ 21 , ..., µ 2b 2 , ..., µ ab a ) . We only implemented balanced designs, i.e., b i = b for all i = 1, . . . , a. Hierarchically nested three-way designs or arbitrary crossed higher-way layouts can be introduced similarly and are implemented as well.
• Repeated Measures and Split Plot Designs are covered by setting d = t, where even hypotheses about within subject factors can be formulated. We exemplify this for profile analyses in the special case of a one-sample RM design with a = 1 as well as for a two-sample RM design with a = 2: Note, that we could also employ more complex factorial structures on the repeated measurements (i.e., more within subject factors) by splitting up the index j.

Examples
To demonstrate the use of the RM and the MANOVA function, we provide several examples for both repeated measures designs and multivariate data in the following. Furthermore, the MANOVA.RM package is equipped with an optional GUI (graphical user interface), based on RGtk2 (Lawrence and Temple Lang, 2010), which will be explained in detail in Section 2.3.3 below.
Both functions are structured similarly: The main input parameters are the formula specifying the outcome on the left hand side and the factor variables of interest on the right, the data and the resampling approach. The latter varies according to the design: the user can choose between a parametric and a wild bootstrap and in the RM design, additionally a permutation approach for the WTS is implemented.

Repeated Measures Designs
The function RM is built as follows: R> RM(formula, data, subject, no.subf = 1, iter = 10000, alpha = 0.05, + resampling = "Perm", CPU, seed, CI.method = "t-quantile", dec = 3) Data need to be provided in long format, i.e., one row per measurement. Here, subject specifies the column name of the subjects variable in the data frame, while no.subf denotes the number of within subject factors considered. Note that in a setting with more than one between subjects factor, the subject ids in the different groups need to be different. Otherwise the program will erroneously assume that these measurements belong to the same subject. The number of cores used for parallel computing as well as a random seed can optionally be specified using CPU and seed, respectively. For calculating the confidence intervals, the user can choose between t-quantiles (the default) and the quantiles based on the resampled WTS. The results are rounded to dec digits.
The function RM returns an object of class "RM" from which the user may obtain plots and summaries of the results using plot(), print() and summary(), respectively. Here, print() returns a short summary of the results, i.e., the values of the test statistics along with degrees of freedom and corresponding p values whereas summary() also displays some descriptive statistics such as the means, sample sizes and confidence intervals for the different factor level combinations. Plotting is based on plotrix (Lemon, 2006). For two-and higher-way layouts, the factors for plotting can be additionally specified in the plot call, see the examples below.

Example 1: One between subject and two within subject factors
For illustration purposes, we consider the data set o2cons, which is included in MANOVA.RM. This data set contains measurements of the oxygen consumption of leukocytes in the presence and absence of inactivated staphylococci at three consecutive time points. Due to the study design, both time and staphylococci are within subject factors while the treatment (Verum vs. Placebo) is a between subject factor (see Friedrich et al., 2017a, for more details).
In this example, we find significant effects of all factors as well as a significant interaction between group and time.

Example 2: Two within subject and two between subject factors
We consider the data set EEG from the MANOVA.RM package: At the Department of Neurology, University Clinic of Salzburg, 160 patients were diagnosed with either Alzheimer's Disease (AD), mild cognitive impairment (MCI), or subjective cognitive complaints without clinically significant deficits (SCC), based on neuropsychological diagnostics (Bathke et al., 2018). This data set contains z-scores for brain rate and Hjorth complexity, each measured at frontal, temporal and central electrode positions and averaged across hemispheres. In addition to standardization, complexity values were multiplied by −1 in order to make them more easily comparable to brain rate values: For brain rate we know that the values decrease with age and pathology, while Hjorth complexity values are known to increase with age and pathology. The three between subject factors considered were sex (men vs. women), diagnosis (AD vs. MCI vs. SCC), and age (< 70 vs. ≥ 70 years). Additionally, the within subject factors region (frontal, temporal, central) and feature (brain rate, complexity) structure the response vector.

Plotting
The RM() function is equipped with a plotting option, displaying the calculated means along with (1 − α) confidence intervals based on t-quantiles. The plot function takes an RM object as argument. In addition, the factor of interest may be specified. If this argument is omitted in a two-or higherway layout, the user is asked to specify the factor for plotting. Furthermore, additional graphical parameters can be used to customize the plots. The optional argument legendpos specifies the position of the legend in higher-way layouts, whereas gap (default 0.1) is the distance introduced between error bars in a higher-way layout. Additionally, the parameter CI.info can be set to TRUE in order to output the means and confidence intervals for the desired interaction.

MANOVA Design
For the analysis of multivariate data, the functions MANOVA and MANOVA.wide are implemented. The difference between the two functions is that the response must be stored in long and wide format for using MANOVA or MANOVA.wide, respectively. The structure of both functions is very similar. They both calculate the WTS for multivariate data in a design with crossed or nested factors. Additionally, the modified ANOVA-type statistic (MATS) is calculated which has the additional advantage of being applicable to designs involving singular covariance matrices and is invariant under scale transformations of the data . The resampling methods provided are a parametric bootstrap approach and a wild bootstrap using Rademacher weights. Note that only balanced nested designs (i.e., the same number of factor levels b for each level of the factor A) with up to three factors are implemented. Designs involving both crossed and nested factors are not implemented. Note that in nested designs, the levels of the nested factor usually have the same labels for all levels of the main factor, i.e., for each level i = 1, . . . , a of the main factor A the nested factor levels are labeled as j = 1, . . . , b i . If the levels of the nested factor are named uniquely, this has to be specified by setting the parameter nested.levels.unique to TRUE.

Data Example MANOVA: Two crossed factors
We again consider the data set EEG from the MANOVA.RM package, but now we ignore the within subject factor structure. Therefore, we are now in a multivariate setting with 6 measurements per patient and three crossed factors sex, age and diagnosis. Due to the small number of subjects in some groups we restrict our analyses to two factors at a time. The analysis of this example is shown below. R> data("EEG") R> EEG_MANOVA <-MANOVA(resp~sex * diagnosis, data = EEG, subject = "id", + resampling = "paramBS", iter = 1000, alpha = 0.01, + CPU = 1, seed = 987) R> summary(EEG_MANOVA) Call: resp~sex * diagnosis  Figure 2: Plot for the interaction between "sex" and "diagnosis" (upper panel) as well as additionally taking "feature" into account (lower panel) in the RM model of the EEG data example.

Wald-Type Statistic (WTS):
Test statistic df p-value sex "12.604" "6" "0.05" diagnosis "55.158" "12" "<0.001" sex:diagnosis "9.79" "12" "0. p-values resampling: paramBS (WTS) paramBS (MATS) sex "0.124" "0.003" diagnosis "<0.001" "<0.001" sex:diagnosis "0.748" "0.223" The output consists of several parts: First, some descriptive statistics of the data set are displayed, namely the sample size and mean for each factor level combination and each dimension (dimensions occur in the same order as in the original data set). In this example, Mean 1 to Mean 3 correspond to the brainrate (temporal, frontal, central) while Mean 4-6 correspond to complexity. Second, the results based on the WTS are displayed. For each factor, the test statistic, degree of freedom and p value is given. For the MATS, only the value of the test statistic is given, since here inference is only based on resampling. The resampling-based p values are finally displayed for both test statistics.
To demonstrate the use of the MANOVA.wide() function, we consider the same data set in wide format, which is also included in the package. In the formula argument, the user now needs to specify the variables of interest bound together via cbind. A subject variable is no longer necessary, as every row of the data set belongs to one patient in wide format data. The output is almost identically to the one obtained from MANOVA with the difference that the mean values are now labeled according to the variable names supplied in the formula argument. R> data("EEGwide") R> EEG_wide <-MANOVA.wide(cbind(brainrate_temporal, brainrate_frontal, + brainrate_central, complexity_temporal, + complexity_frontal, complexity_central)~sex * diagnosis, + data = EEGwide, resampling = "paramBS", iter = 1000, + alpha = 0.01, CPU = 1, seed = 987) R> summary(EEG_wide) Call: cbind(brainrate_temporal, brainrate_frontal, brainrate_central, complexity_temporal, complexity_frontal, complexity_central) p-values resampling: paramBS (WTS) paramBS (MATS) sex "0.122" "0.005" diagnosis "<0.001" "<0.001" sex:diagnosis "0.742" "0.21" In this example, MATS detects a significant effect of sex, a finding that is not shared by the p value based on the parametric bootstrap WTS.

Confidence Regions
The MANOVA functions are equipped with a function for calculating and plotting of confidence regions. Details on the methods can be found in . We would like to point out, however, that the MATS-based confidence regions have to be interpreted differently from the more usual WTS-based ones. For the latter, the WTS is compared to a fixed critical value (in the asymptotic choice from a χ 2 -distribution) and we thus obtain geometric ellipsoids as the WTS is more or less a Mahalanobis distance in the inverse covariance matrix. The MATS statistic only involves the variances of the covariance matrix and we thus obtain a different geometric shape of the corresponding confidence region. However, here the correlation is implicitly involved in the critical value which now (different to the WTS) depends on the covariance matrix. More precisely, a confidence region for the vector of contrasts Hµ based on the MATS is determined by the set of all Hµ such that where c * 1−α denotes the quantile of the respective resampling distribution. A confidence ellipsoid is now obtained based on the eigenvalues λ s and eigenvectors e s of H D N H . The classical WTS-based confidence ellipsoids, in contrast, are based on eigenvectors and eigenvalues of H Σ N H instead. Confidence regions can be calculated using the conf.reg function. Note that confidence regions can only be plotted in designs with 2 dimensions.

R> conf.reg(object, nullhypo)
Object must be an object of class "MANOVA", i.e., created using either MANOVA or MANOVA.wide, whereas nullhypo specifies the desired null hypothesis, i.e., the contrast of interest in designs involving more than one factor. As an example, we consider the data set water from the HSAUR package (Everitt and Hothorn, 2017). The data set contains measurements of mortality and drinking water hardness for 61 cities in England and Wales. Suppose we want to analyse whether these measurements differ between northern and southern towns. Since the data set is in wide format, we need to use the MANOVA.wide function. R> library("HSAUR") R> data("water") R> test <-MANOVA.wide(cbind(mortality, hardness)~location, data = water, + iter = 1000, resampling = "paramBS", CPU = 1, seed = 123) R> summary(test) R> cr <-conf.reg(test) R> cr R> plot(cr, xlab = "Difference in mortality", + ylab ="Difference in water hardness") We find significant differences in mortality and water hardness between northern and southern towns.
The confidence region is returned as an ellipsoid specified by its center as well as its axes, which extend Scale units into the direction of the respective eigenvector. For two-dimensional outcomes as in this example, the confidence ellipsoid can also be plotted using the ellipse package (Murdoch and Chow, 2018), see Figure 3.

Nested design
To create a data example for a nested design, we use the curdies data set from the GFD package and extend it by introducing an artificial second outcome variable. In this data set, the levels of the nested factor (site) are named uniquely, i.e., levels 1-3 of factor site belong to "WINTER", whereas levels 4-6 belong to "SUMMER". Therefore, nested.levels.unique must be set to TRUE. The code for the analysis using both wide and long format is presented below. Calculation of simultaneous confidence intervals and p values for contrasts of the mean vector is based on the sum statistic, see  for details. Confidence intervals are calculated based on summary effects, i.e., averaging over all dimensions, whereas the returned pvalues are multivariate. Note that the confidence intervals and p values returned are simultaneous, i. e., they maintain the given alpha-level. Such contrasts include, e. g., Tukey's all-pairwise comparisons or Dunnett's many-to-one comparisons, see e. g. Hothorn et al. (2008a) for more examples. Confidence intervals for contrasts of the mean vector can be calculated using the function simCI, which is build on contrMat from the multcomp package (Hothorn et al., 2008b): simCI(object, contrast, contmat, type, base) Here, object is an object of class "MANOVA". The user can choose between pairwise or user-defined contrasts. For user-defined constrast (contrast = "user-defined"), the contrast matrix of interest must be specified in contmat. Pairwise comparisons (contrast = "pairwise") are calculated using the contrMat function of multcomp and accordingly, type and base specify the type of the pairwise comparison and the baseline group for Dunnett contrasts, see Hothorn et al. (2008b) for details on these parameters. To exemplify its application we reconsider the EEG example from above: R> # pairwise comparison using Tukey contrasts R> simCI(EEG_MANOVA, contrast = "pairwise", type = "Tukey")  The output is two-fold: First, the multivariate p values for the desired contrasts are returned. The second part of the output provides simultaneous confidence intervals for summary effects by averaging over all dimensions.
If the global null hypothesis, e. g.
has been rejected, it is usually of interest to further investigate which univariate endpoints caused the rejection. A straight-forward way to do this is to calculate the univariate p values and adjust them for multiple testing, e. g., using Bonferroni correction. Consider the one-way layout above: Since the global null hypothesis can be rejected, we now wish to analyze which of the two univariate endpoints caused this rejection. R> EEG1 <-MANOVA.wide(brainrate_temporal~diagnosis, data = EEGwide, + iter = 1000, seed = 987, CPU = 1) R> EEG2 <-MANOVA.wide(brainrate_central~diagnosis, data = EEGwide, + iter = 1000, seed = 987, CPU = 1) R> p.adjust(c(EEG1$resampling[, 2], EEG2$resampling[, 2]), + method = "bonferroni") [1] 0 0 Thus, in this example both endpoints showed significant effects. Note that it is often possible to conduct post-hoc comparisons according to the closure principle, thus avoiding the need to correct for multiple comparisons. Implementation of these methods for both MANOVA and RM designs is part of future research.

Graphical user interface
The GUI is started in R with the command GUI.RM(), GUI.MANOVA() and GUI.MANOVAwide() for repeated measures designs and multivariate data in long or wide format, respectively. Note that the GUI depends on RGtk2 and will only work if RGtk2 is installed. The user can specify the data location (either directly or via the "load data" button) and the formula as well as the number of iterations, the significance level α, the number of within subject factors (for repeated measures designs) and the name of the subject variable, see Figure 4. Furthermore, the user has the choice between the three resampling approaches "Perm" (only for RM designs), "paramBS" and "WildBS" denoting the permutation procedure, the parametric bootstrap and the wild bootstrap, respectively. Additionally, one can specify whether or not headers are included in the data file, and which separator and character symbols are used for decimals in the data file. The GUI for repeated measures also provides a plotting option, which generates a new window for specifying the factors to be plotted (in higher-way layouts) along with a few plotting parameters, see Figure 5.   The left GUI is for the one-way layout (no choice of factors possible), the right one is for a two-way layout with an example for plotting interactions.

Discussion and Outlook
We have explicitly described the usage of the R package MANOVA.RM for analyzing various nonparametric multivariate MANOVA and RM designs making use of novel bootstrap-and permutationapproaches. Moreover, the corresponding models and inference procedures that have been derived and theoretically analyzed in previous papers are explained as well. In particular, three different test statistics of Wald-, ANOVA-and modified ANOVA-type are implemented together with appropriate critical values derived from asymptotic considerations, approximations or novel resampling approaches. Here, the latter is recommended in case of small to moderate sample sizes. All methods can be applied without assuming usual presumptions such as multivariate normality or specific covariance structures. Moreover, all procedures are particularly constructed to tackle covariance matrix heterogeneity across groups or even covariance singularity (in case of the MATS). In this way MANOVA.RM provides a flexible tool box for inferring hypotheses about (i) main and interaction effects in general factorial MANOVA and (ii) between and within subject effects in RM designs with possibly complex factorial structures on both, between and within subject factors.
In addition, we have placed a graphical user interface (GUI) at the users disposal to allow for a simple and intuitive use. It is planned to update the package on a regular basis; respecting the development of new procedures for general RM and MANOVA designs. For example, our working group is currently investigating the implementation of covariates in the above model in theoretical research and the resulting procedure may be incorporated in the future. Other topics include the possible implementation/improvement of subsequent multiple comparisons by the closure principle, as in Burchett et al. (2017), for both MANOVA and RM designs.