1 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:
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.

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 subplot factors (e.g. time) are of interest. On the other hand, MANOVA settings are usually only designed to detect effects of the observed factors (and interactions thereof) on the multivariate outcome vectors.
Despite their differences, MANOVA and RMtype techniques share the same advantages over classical univariate endpointwise—ANOVAtype—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.
Focusing on metric data and meanbased procedures, MANOVA and RMmodels are typically inferred by means of “classical” procedures such as Wilks’ Lambda, LawleyHotelling, Roy’s largest root (Davis, 2002; Johnson and Wichern, 2007; Anderson, 2001)
or (generalized) linear mixed models with generalized estimating equations. For the classical oneway layout, these methods are implemented in
R within the manova function in the stats package, where one can choose between the options ‘Pillai‘, ‘Wilks‘, ‘HotellingLawley‘ and ‘Roy‘. Nonparametric rankbased methods for null hypotheses formulated in distribution functions are implemented within the packages npmv for one and twoway 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 Waldtype test in the univariate setting, can also be used.(Generalized) linear mixed models are implemented in the lm and the glm function (package stats, R Core Team, 2016) for univariate data as well as in the SCGLR package for Generalized Linear Model estimation in the context of multivariate data (Cornu et al., 2016). The Anova and Manova function in the car package (Fox and Weisberg, 2011) calculate typeII and typeIII analysisofvariance 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 typeIerrors, cf.
Vallejo et al. (2001); Lix and Keselman (2004); Vallejo Seco et al. (2007); LivacicRojas 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); Zhang (2011); Harrar and Bathke (2012); Konietschke et al. (2015); Xiao and Zhang (2016); Bathke et al. (2016); McFarquhar et al. (2016); Friedrich et al. (2017a); LivacicRojas et al. (2017); Friedrich and Pauly (2017)and the references cited therein. Here, we focus on statistical methods that are valid in the multivariate BehrensFisher situation—equal covariance matrices across the groups is not assumed—and provide accurate inferential results in terms of pvalue estimates and confidence intervals for the parameters of interest. In particular, general Waldtype test statistics (for MANOVA and RM), ANOVAtype statistics (for RM) and modified ANOVAtype tests (for MANOVA) are implemented in
MANOVA.RM 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 are appropriate methods in the BehrensFisher situation.
To make the methods freely accessible we have provided the R package MANOVA.RM for daily statistical analyses. It is available from the R Archive at
https://CRAN.Rproject.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 (R package stats, R Core Team, 2016). Its userfriendly application not only provides the values and test statistics of interest but also a descriptive overview together with componentwise twosided confidence intervals. Moreover, the MANOVA function even allows for an easy calculation and confidence ellipsoids plots for specified multivariate contrasts as described in Friedrich and Pauly (2017).
Specifically, for testing multivariate main and interaction effects in one, two and higherway MANOVA models, the MANOVA function provides the

Waldtype statistic (WTS) proposed by Konietschke et al. (2015) using a parametric bootstrap, a wild bootstrap or its asymptotic distribution for value computations, and

the modified ANOVAtype statistic (MATS) proposed by Friedrich and Pauly (2017) using a parametric or wild bootstrap procedure for value computations.
In addition to multivariate groupwise effects, the RM function also allows to test hypotheses formulated across subplot factors. The implemented test statistics are

the ANOVAtype statistic (ATS) using an approximation as considered in Brunner (2001) as well as a parametric and a wild bootstrap approach and
The paper is organized as follows: In Section 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 3. Finally, the paper closes with a discussion in Section 4.
Throughout the paper we use the subsequent notation from multivariate linear models: For we denote by the dimensional centering matrix, by the
dimensional identity matrix and by
the matrix of 1’s, i.e., , where is the dimensional column vector of 1’s.2 Statistical model and inference methods
For both the RM and the MANOVA design being equipped with an arbitrary number of fixed factors, we consider the general linear model given by variate random vectors
(1) 
Here, denotes the experimental unit or subject in group . Note, that a higherway factorial structure on the groups/wholeplots or subplots can be achieved by subindexing the indices (group/wholeplot) or (subplot) into or . In this model is the mean vector in group and for each fixed it is assumed that the error terms are independent and identically distributed variate random vectors with mean and existing variances . For the WTStype procedures we additionally assume positive definite covariance matrices
and existing finite fourth moments
In this model, hypotheses for RM or MANOVA can be formulated by means of an adequate contrast hypothesis matrix by
Let denote the vector of pooled group means and the estimated covariance of . Here, and . In this setup Konietschke et al. (2015) propose
a Waldtype statistic (WTS)
(2) 
for testing , where , , and denotes the MoorePenrose inverse of the matrix . Since its asymptotic null distribution provides a poor finite sample approximation, they propose the following asymptotic modelbased bootstrap approach: Given the data let be independent random vectors that are used for recalculating the test statistic as , where . Denoting by the corresponding
quantile of the (conditional) distribution of
the test rejects if . 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. (2016)). However, Friedrich et al. (2017a) proposed a more favourable technique for RM designs. 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 their permutation test for RM models rejects if . Here denotes the quantile of the (conditional) distribution of the permutation version of the test statistic . As shown in extensive simulations in Friedrich et al. (2017a) and the corresponding supporting information this ’wild’ permutation WTS method controls the type1 error rate very well. Note that this procedure is only applicable for RM due to the commensurate nature of their components. In MANOVA setups 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 wellestablished ANOVAtype statistic (ATS)
(3) 
by Brunner (2001) is implemented together with the enhanced approximation of the statistic proposed by Brunner et al. (1997, 2002) and also 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 or since it waives the MoorePenrose inverse involved in (2).
Similar to the permuted WTS the ATS given in (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 or , Friedrich and Pauly (2017) have recently proposed the novel MATS (modified ATS)
(4) 
Here, the involved diagonal matrix of the empirical variances of component in group , deduces an invariance under componentwise scale transformations of the MATS. To obtain an accurate finite sample performance, it is also equipped with an asymptotic model based bootstrap approach. That is, MATS rejects if , where is the quantile of the (conditional) distribution of the bootstrapped statistic . In addition, we implemented a wild bootstrap approach, which is based on multiplying the centered data vectors with random weights fulfilling and . In the package, we implemented Rademacher distributed weights, i.e., random signs. Extensive simulations in Friedrich and Pauly (2017) 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.
2.1 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 being implemented in MANOVA.RM are even applicable in higherway layouts than being presented here; and the list should not be seen as the limited application of the package. The models are derived by subindexing the index in model (1) in the following ways:

OneWay (): Writing we have with
and obtain the null hypothesis of ’no group’ or ’factor
’ effect asIn case of this includes the famous multivariate BehrensFisher 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 TwoWay (): Splitting the index into two and writing we obtain the model with The corresponding null hypotheses of no main effects in or and no interaction effect between and can be written as:

Hierarchically nested TwoWay (): A fixed subcategory within factor can be introduced via the model with Here, the hypotheses of no main effect or no subcategory main effect can be written as
with and .
We only implemented balanced designs, i.e., for all . Hierarchically nested threeway designs or arbitrary crossed higherway layouts can be introduced similarly and are implemented as well.

Repeated Measures and Split Plot Designs are covered by setting , where even hypotheses about subplots can be formulated. We exemplify this for profile analyses in the special case of a onesample RM design with
as well as for a twosample RM design with :
with , and .
parallel flat identical
Note, that we could also employ more complex factorial structures on the repeated measurements (i.e., more subplot factors) by splitting up the index .
3 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, based on RGtk2 (Lawrence and Temple Lang, 2010), which will be explained in detail in Section 3.3 below.
3.1 Repeated Measures Designs
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
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 higherway layouts, the factors for plotting can be additionally specified in the plot call, see the examples below. R¿ RM(formula, data, subject, no.subf = 1, iter = 10000, + alpha = 0.05, resampling = ”Perm”, CPU, seed, + CI.method = ”tquantile”, 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 subplot factors considered. 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 quantiles (the default) and the quantiles based on the resampled WTS. The results are rounded to dec digits.
3.1.1 Example 1: One wholeplot and two subplot 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 subplot factors while the treatment (Verum vs. Placebo) is a wholeplot factor (Friedrich et al., 2017a).
R¿ data(”o2cons”) R¿ model1 ¡ RM(O2 Group * Staphylococci * Time, data = o2cons, + subject = ”Subject”, no.subf = 2, iter = 1000, + resampling = ”Perm”, seed = 1234) R¿ summary(model1)
Call: O2 Group * Staphylococci * Time
Descriptive: Group Staphylococci Time n Means Lower 95 1 P 0 6 12 1.322 1.150 1.493 5 P 0 12 12 2.430 2.196 2.664 9 P 0 18 12 3.425 3.123 3.727 3 P 1 6 12 1.618 1.479 1.758 7 P 1 12 12 2.434 2.164 2.704 11 P 1 18 12 3.527 3.273 3.781 2 V 0 6 12 1.394 1.201 1.588 6 V 0 12 12 2.570 2.355 2.785 10 V 0 18 12 3.677 3.374 3.979 4 V 1 6 12 1.656 1.471 1.840 8 V 1 12 12 2.799 2.500 3.098 12 V 1 18 12 4.029 3.802 4.257
WaldType Statistic (WTS): Test statistic df pvalue Group 11.167 1 0.001 Staphylococci 20.401 1 0.000 Group:Staphylococci 2.554 1 0.110 Time 4113.057 2 0.000 Group:Time 24.105 2 0.000 Staphylococci:Time 4.334 2 0.115 Group:Staphylococci:Time 4.303 2 0.116
ANOVAType Statistic (ATS): Test statistic df1 df2 pvalue Group 11.167 1.000 316.278 0.001 Staphylococci 20.401 1.000 Inf 0.000 Group:Staphylococci 2.554 1.000 Inf 0.110 Time 960.208 1.524 Inf 0.000 Group:Time 5.393 1.524 Inf 0.009 Staphylococci:Time 2.366 1.983 Inf 0.094 Group:Staphylococci:Time 2.147 1.983 Inf 0.117
pvalues resampling: Perm (WTS) Perm (ATS) Group 0.005 NA Staphylococci 0.000 NA Group:Staphylococci 0.138 NA Time 0.000 NA Group:Time 0.000 NA Staphylococci:Time 0.161 NA Group:Staphylococci:Time 0.163 NA
The output consists of four parts: model1$Descriptive gives an overview of the descriptive statistics: The number of observations, mean and confidence intervals are displayed for each factor level combination. Second, model1$WTS contains the results for the Waldtype test: The test statistic, degree of freedom and values based on the asymptotic distribution are displayed. Note that the approximation is very liberal for small sample sizes, cf. Konietschke et al. (2015); Friedrich et al. (2017a). The corresponding results based on the ATS are contained within model1$ATS. This test statistic tends to rather conservative decisions in case of small sample sizes and is even asymptotically only an approximation, thus not providing an asymptotic level test, see Brunner (2001); Friedrich et al. (2017a). Finally, model1$resampling contains the values based on the chosen resampling approach. For the ATS, the permutation approach is not feasible since it would result in an incorrect covariance structure, and is therefore not implemented. Due to the above mentioned issues for small sample sizes, the respective resampling procedure is recommended for such situations.
In this example, we find significant effects of all factors as well as a significant interaction between group and time.
3.1.2 Example 2: Two subplot and two wholeplot 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 AD, MCI, or SCC, based on neuropsychological diagnostics (Bathke et al., 2016). This data set contains 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 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 wholeplot factors considered were sex (men vs. women), diagnosis (AD vs. MCI vs. SCC), and age ( vs. years). Additionally, the subplot factors region (frontal, temporal, central) and feature (brain rate, complexity) structure the response vector.
R¿ data(”EEG”) R¿ EEG_model ¡ RM(resp sex * diagnosis * feature * region, + data = EEG, subject = ”id”, no.subf = 2, + resampling = ”WildBS”, iter = 1000, alpha = 0.01, + CPU = 4, seed = 123) R¿ summary(EEG_model)
Call: resp sex * diagnosis * feature * region
Descriptive: sex diagnosis feature region n Means Lower 99 1 M AD brainrate central 12 1.010 4.881 2.861 13 M AD brainrate frontal 12 1.007 4.991 2.977 25 M AD brainrate temporal 12 0.987 4.493 2.519 7 M AD complexity central 12 1.488 10.053 7.077 19 M AD complexity frontal 12 1.086 6.906 4.735 31 M AD complexity temporal 12 1.320 7.203 4.562 3 M MCI brainrate central 27 0.447 1.591 0.696 15 M MCI brainrate frontal 27 0.464 1.646 0.719 27 M MCI brainrate temporal 27 0.506 1.584 0.572 9 M MCI complexity central 27 0.257 1.139 0.625 21 M MCI complexity frontal 27 0.459 1.997 1.079 33 M MCI complexity temporal 27 0.490 1.796 0.816 5 M SCC brainrate central 20 0.459 0.414 1.332 17 M SCC brainrate frontal 20 0.243 0.670 1.156 29 M SCC brainrate temporal 20 0.409 1.210 2.028 11 M SCC complexity central 20 0.349 0.070 0.767 23 M SCC complexity frontal 20 0.095 1.037 1.227 35 M SCC complexity temporal 20 0.314 0.598 1.226 2 W AD brainrate central 24 0.294 1.978 1.391 14 W AD brainrate frontal 24 0.159 1.813 1.495 26 W AD brainrate temporal 24 0.285 1.776 1.206 8 W AD complexity central 24 0.128 1.372 1.116 20 W AD complexity frontal 24 0.026 1.212 1.264 32 W AD complexity temporal 24 0.194 1.670 1.283 4 W MCI brainrate central 30 0.106 1.076 0.863 16 W MCI brainrate frontal 30 0.074 1.032 0.885 28 W MCI brainrate temporal 30 0.069 1.064 0.925 10 W MCI complexity central 30 0.094 0.464 0.652 22 W MCI complexity frontal 30 0.131 0.768 1.031 34 W MCI complexity temporal 30 0.121 0.652 0.895 6 W SCC brainrate central 47 0.537 0.049 1.124 18 W SCC brainrate frontal 47 0.548 0.062 1.159 30 W SCC brainrate temporal 47 0.559 0.015 1.133 12 W SCC complexity central 47 0.384 0.110 0.659 24 W SCC complexity frontal 47 0.403 0.038 0.845 36 W SCC complexity temporal 47 0.506 0.132 0.880
WaldType Statistic (WTS): Test statistic df pvalue sex 9.973 1 0.002 diagnosis 42.383 2 0.000 sex:diagnosis 3.777 2 0.151 feature 0.086 1 0.769 sex:feature 2.167 1 0.141 diagnosis:feature 5.317 2 0.070 sex:diagnosis:feature 1.735 2 0.420 region 0.070 2 0.966 sex:region 0.876 2 0.645 diagnosis:region 6.121 4 0.190 sex:diagnosis:region 1.532 4 0.821 feature:region 0.652 2 0.722 sex:feature:region 0.423 2 0.810 diagnosis:feature:region 7.142 4 0.129 sex:diagnosis:feature:region 2.274 4 0.686
ANOVAType Statistic (ATS): Test statistic df1 df2 pvalue sex 9.973 1.000 657.416 0.002 diagnosis 13.124 1.343 657.416 0.000 sex:diagnosis 1.904 1.343 657.416 0.164 feature 0.086 1.000 Inf 0.769 sex:feature 2.167 1.000 Inf 0.141 diagnosis:feature 1.437 1.562 Inf 0.238 sex:diagnosis:feature 1.031 1.562 Inf 0.342 region 0.018 1.611 Inf 0.965 sex:region 0.371 1.611 Inf 0.644 diagnosis:region 1.091 2.046 Inf 0.337 sex:diagnosis:region 0.376 2.046 Inf 0.691 feature:region 0.126 1.421 Inf 0.810 sex:feature:region 0.077 1.421 Inf 0.864 diagnosis:feature:region 0.829 1.624 Inf 0.415 sex:diagnosis:feature:region 0.611 1.624 Inf 0.510
pvalues resampling: WildBS (WTS) WildBS (ATS) sex 0.000 0.000 diagnosis 0.000 0.000 sex:diagnosis 0.119 0.124 feature 0.798 0.798 sex:feature 0.152 0.152 diagnosis:feature 0.067 0.249 sex:diagnosis:feature 0.445 0.362 region 0.967 0.980 sex:region 0.691 0.728 diagnosis:region 0.182 0.338 sex:diagnosis:region 0.863 0.814 feature:region 0.814 0.926 sex:feature:region 0.881 0.951 diagnosis:feature:region 0.098 0.519 sex:diagnosis:feature:region 0.764 0.683
We find significant effects at level of the wholeplot factors sex and diagnosis, while none of the subplot factors or interactions become significant.
3.1.3 Plotting
The RM() function is equipped with a plotting option, displaying the calculated means along with confidence intervals based on 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 higherway layouts, whereas gap (default 0.1) is the distance introduced between error bars in a higherway layout. R¿ plot(EEG_model, factor = ”sex”, main = + ”Effect of sex on EEG values”) R¿ plot(EEG_model, factor = ”sex:diagnosis”, legendpos = ”topleft”, + col = c(4, 2), ylim = c(1.8, 0.8)) R¿ plot(EEG_model, factor = ”sex:diagnosis:feature”, + legendpos = ”bottomright”, gap = 0.05)
3.2 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 ANOVAtype 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 (Friedrich and Pauly, 2017). 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 for each level of the factor ) 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 of the main factor the nested factor levels are labeled as . If the levels of the nested factor are named uniquely, this has to be specified by setting the parameter nested.levels.unique to TRUE.
R¿ MANOVA(formula, data, subject, iter = 10000, alpha = 0.05, resampling = ”paramBS”, CPU, seed, nested.levels.unique = FALSE, dec = 3) R¿ MANOVA.wide(formula, data, iter = 10000, alpha = 0.05, resampling = ”paramBS”, CPU, seed, nested.levels.unique = FALSE, dec = 3)
The only difference between MANOVA and MANOVA.wide in the function call except from the different shape of the formula (see examples below) is the subject variable, which needs to be specified for MANOVA only.
3.2.1 Data Example Manova: Two crossed factors
We again consider the data set EEG from the MANOVA.RM package, but now we ignore the subplot 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 (e.g., only 2 male patients aged 70 were diagnosed with AD) 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
Descriptive: sex diagnosis n Mean 1 Mean 2 Mean 3 Mean 4 Mean 5 Mean 6 1 M AD 12 0.987 1.007 1.010 1.320 1.086 1.488 3 M MCI 27 0.506 0.464 0.447 0.490 0.459 0.257 5 M SCC 20 0.409 0.243 0.459 0.314 0.095 0.349 2 W AD 24 0.285 0.159 0.294 0.194 0.026 0.128 4 W MCI 30 0.069 0.074 0.106 0.121 0.131 0.094 6 W SCC 47 0.559 0.548 0.537 0.506 0.403 0.384
WaldType Statistic (WTS): Test statistic df pvalue sex 12.604 6 0.050 diagnosis 55.158 12 0.000 sex:diagnosis 9.790 12 0.634
modified ANOVAType Statistic (MATS): Test statistic sex 45.263 diagnosis 194.165 sex:diagnosis 18.401
pvalues resampling: paramBS (WTS) paramBS (MATS) sex 0.124 0.003 diagnosis 0.000 0.000 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 value is given. For the MATS, only the value of the test statistic is given, since here inference is only based on resampling. The resamplingbased 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) sex * diagnosis
Descriptive: sex diagnosis n brainrate_temporal brainrate_frontal brainrate_central 1 M AD 12 0.987 1.007 1.010 2 W AD 27 0.506 0.464 0.447 3 M MCI 20 0.409 0.243 0.459 4 W MCI 24 0.285 0.159 0.294 5 M SCC 30 0.069 0.074 0.106 6 W SCC 47 0.559 0.548 0.537 complexity_temporal complexity_frontal complexity_central 1 1.320 1.086 1.488 2 0.490 0.459 0.257 3 0.314 0.095 0.349 4 0.194 0.026 0.128 5 0.121 0.131 0.094 6 0.506 0.403 0.384
WaldType Statistic (WTS): Test statistic df pvalue sex 12.604 6 0.050 diagnosis 55.158 12 0.000 sex:diagnosis 9.790 12 0.634
modified ANOVAType Statistic (MATS): Test statistic sex 45.263 diagnosis 194.165 sex:diagnosis 18.401
pvalues resampling: paramBS (WTS) paramBS (MATS) sex 0.124 0.003 diagnosis 0.000 0.000 sex:diagnosis 0.748 0.223
In this example, MATS detects a significant effect of sex, a finding that is not shared by the value based on the parametric bootstrap WTS.
3.2.2 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 Friedrich and Pauly (2017). 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, 2015). 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)
Call: cbind(mortality, hardness) location
Descriptive: location n mortality hardness North North 35 1633.600 30.400 South South 26 1376.808 69.769
WaldType Statistic (WTS): Test statistic df pvalue 51.584 2.000 0.000
modified ANOVAType Statistic (MATS): [,1] [1,] 69.882
pvalues resampling: paramBS (WTS) paramBS (MATS) 0 0
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 twodimensional outcomes as in this example, the confidence ellipsoid can also be plotted, see Figure
3.Center: [,1] [1,] 256.792 [2,] 39.369
Scale: [1] 10.852716 2.736354
Eigenvectors: [,1] [,2] [1,] 1 0 [2,] 0 1
3.2.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 13 of factor site belong to ”WINTER”, whereas levels 46 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.
R¿ library(GFD) R¿ data(curdies) R¿ set.seed(123) R¿ curdiesdugesia + rnorm(36)
R¿ # first possibility: MANOVA.wide R¿ fit1 ¡ MANOVA.wide(cbind(dugesia, dug2) season + season:site, + data = curdies, iter = 100, nested.levels.unique = TRUE, + seed = 123, CPU = 1)
R¿ # second possibility: MANOVA (long format) R¿ dug ¡ c(curdiesdug2) R¿ season ¡ rep(curdiessite, 2) R¿ curd ¡ data.frame(dug, season, site, subject = rep(1:36, 2)) R¿ fit2 ¡ MANOVA(dug season + season:site, data = curd, + subject = ”subject”, nested.levels.unique = TRUE, + seed = 123, iter = 100, CPU = 1)
R¿ # comparison of results R¿ summary(fit1) R¿ summary(fit2)
Call: cbind(dugesia, dug2) season + season:site
Descriptive: season site n dugesia dug2 1 SUMMER 4 6 0.419 0.050 2 SUMMER 5 6 0.229 0.028 3 SUMMER 6 6 0.194 0.763 4 WINTER 1 6 2.049 2.497 5 WINTER 2 6 4.182 4.123 6 WINTER 3 6 0.678 0.724
WaldType Statistic (WTS): Test statistic df pvalue season 6.999 2 0.030 season:site 16.621 8 0.034
modified ANOVAType Statistic (MATS): Test statistic season 12.296 season:site 15.064
pvalues resampling: paramBS (WTS) paramBS (MATS) season 0.04 0.04 season:site 0.28 0.18
Call: dug season + season:site
Descriptive: season site n Mean 1 Mean 2 1 SUMMER 4 6 0.419 0.050 2 SUMMER 5 6 0.229 0.028 3 SUMMER 6 6 0.194 0.763 4 WINTER 1 6 2.049 2.497 5 WINTER 2 6 4.182 4.123 6 WINTER 3 6 0.678 0.724
WaldType Statistic (WTS): Test statistic df pvalue season 6.999 2 0.030 season:site 16.621 8 0.034
modified ANOVAType Statistic (MATS): Test statistic season 12.296 season:site 15.064
pvalues resampling: paramBS (WTS) paramBS (MATS) season 0.04 0.04 season:site 0.28 0.18
3.3 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 subplot 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 higherway layouts) along with a few plotting parameters, see Figure 5.
R¿ library(”MANOVA.RM”) R¿ GUI.RM()
4 Discussion and Outlook
We have explicitly described the usage of the R package MANOVA.RM for analyzing various multivariate MANOVA and RM designs. 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 ANOVAtype are implemented together with appropriate critical values derived from asymptotic considerations, approximations or resampling. 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) whole and subplot effects in RM designs with possibly complex factorial structures on both, whole and subplots.
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 of subsequent multiple comparisons, e.g. by the closure principle.
Acknowledgments
The work of Sarah Friedrich and Markus Pauly was supported by the German Research Foundation project DFGPA 2409/31.
References
 Anderson (2001) Anderson, M. J. (2001). A new method for nonparametric multivariate analysis of variance. Austral ecology, 26(1):32–46.
 Bathke et al. (2016) Bathke, A., Friedrich, S., Konietschke, F., Pauly, M., Staffen, W., Strobl, N., and Höller, Y. (2016). Using EEG, SPECT, and multivariate resampling methods to differentiate between Alzheimer’s and other cognitive impairments. arXiv preprint arXiv:1606.09004.
 Brunner (2001) Brunner, E. (2001). Asymptotic and approximate analysis of repeated measures designs under heteroscedasticity. Mathematical Statistics with Applications in Biometry.
 Brunner et al. (1997) Brunner, E., Dette, H., and Munk, A. (1997). Boxtype approximations in nonparametric factorial designs. Journal of the American Statistical Association, 92(440):1494–1502.
 Brunner et al. (2002) Brunner, E., Domhof, S., and Langer, F. (2002). Nonparametric Analysis of Longitudinal Data in Factorial Experiments. Wiley, New York, USA.
 Burchett et al. (2017) Burchett, W. W., Ellis, A. R., Harrar, S. W., and Bathke, A. C. (2017). Nonparametric inference for multivariate data: The R package npmv. Journal of Statistical Software, 76(4):1–18.
 Christensen and Rencher (1997) Christensen, W. F. and Rencher, A. C. (1997). A comparison of Type I error rates and power levels for seven solutions to the multivariate BehrensFisher problem. Communications in Statistics  Simulation and Computation, 26(4):1251–1273.

Cornu et al. (2016)
Cornu, G., Mortier, F., Trottier, C., and Bry, X. (2016).
SCGLR: Supervised Component Generalized Linear Regression
. R package version 2.0.3.  Davis (2002) Davis, C. S. (2002). Statistical methods for the analysis of repeated measurements. Springer, New York.
 Everitt and Hothorn (2015) Everitt, B. S. and Hothorn, T. (2015). HSAUR: A Handbook of Statistical Analyses Using R (1st Edition). R package version 1.37.
 Fox and Weisberg (2011) Fox, J. and Weisberg, S. (2011). An R Companion to Applied Regression. Sage, Thousand Oaks CA, second edition.
 Friedrich et al. (2017a) Friedrich, S., Brunner, E., and Pauly, M. (2017a). Permuting longitudinal data in spite of the dependencies. Journal of Multivariate Analysis, 153:255–265.
 Friedrich et al. (2017b) Friedrich, S., Konietschke, F., and Pauly, M. (2017b). GFD: An R package for the analysis of general factorial designs. Journal of Statistical Software, Code Snippets, 79(1):1–18.
 Friedrich and Pauly (2017) Friedrich, S. and Pauly, M. (2017). MATS: Inference for potentially singular and heteroscedastic MANOVA. arXiv preprint arXiv:1704.03731.
 Gupta et al. (2008) Gupta, A. K., Harrar, S. W., and Fujikoshi, Y. (2008). MANOVA for large hypothesis degrees of freedom under nonnormality. Test, 17(1):120–137.
 Harrar and Bathke (2012) Harrar, S. W. and Bathke, A. C. (2012). A modified twofactor multivariate analysis of variance: asymptotics and small sample approximations. Annals of the Institute of Statistical Mathematics, 64(1):135–165.
 Johnson and Wichern (2007) Johnson, R. A. and Wichern, D. W. (2007). Applied multivariate statistical analysis. Prentice Hall.
 Konietschke et al. (2015) Konietschke, F., Bathke, A., Harrar, S., and Pauly, M. (2015). Parametric and nonparametric bootstrap methods for general MANOVA. Journal of Multivariate Analysis, 140:291–301.

Krishnamoorthy and Yu (2004)
Krishnamoorthy, K. and Yu, J. (2004).
Modified Nel and Van der Merwe test for the multivariate
Behrens–Fisher problem.
Statistics & probability letters
, 66(2):161–169.  Lawrence and Temple Lang (2010) Lawrence, M. and Temple Lang, D. (2010). RGtk2: A graphical user interface toolkit for R. Journal of Statistical Software, 37(8):1–52.
 Lemon (2006) Lemon, J. (2006). Plotrix: A package in the red light district of R. RNews, 6(4):8–12.
 LivacicRojas et al. (2010) LivacicRojas, P., Vallejo, G., and Fernandez, P. (2010). Analysis of type i error rates of univariate and multivariate procedures in repeated measures designs. Communications in Statistics  Simulation and Computation, 39(3):624–640.
 LivacicRojas et al. (2017) LivacicRojas, P., Vallejo, G., Fernández, P., and TueroHerrero, E. (2017). Power of modified brownforsythe and mixedmodel approaches in splitplot designs. Methodology.
 Lix and Keselman (2004) Lix, L. M. and Keselman, H. (2004). Multivariate tests of means in independent groups designs: Effects of covariance heterogeneity and nonnormality. Evaluation & the health professions, 27(1):45–69.
 Lix and Lloyd (2007) Lix, L. M. and Lloyd, A. M. (2007). A comparison of procedures for the analysis of multivariate repeated measurements. Journal of Modern Applied Statistical Methods, 6(2):380–398.
 McFarquhar et al. (2016) McFarquhar, M., McKie, S., Emsley, R., Suckling, J., Elliott, R., and Williams, S. (2016). Multivariate and repeated measures (mrm): A new toolbox for dependent and multimodal grouplevel neuroimaging data. NeuroImage, 132:373–389.
 Nel and Van der Merwe (1986) Nel, D. and Van der Merwe, C. (1986). A solution to the multivariate BehrensFisher problem. Communications in StatisticsTheory and Methods, 15(12):3719–3735.
 Noguchi et al. (2012) Noguchi, K., Gel, Y. R., Brunner, E., and Konietschke, F. (2012). nparLD: An R software package for the nonparametric analysis of longitudinal data in factorial experiments. Journal of Statistical Software, 50(12):1–23.
 R Core Team (2016) R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
 Vallejo et al. (2001) Vallejo, G., Fidalgo, A., and Fernandez, P. (2001). Effects of covariance heterogeneity on three procedures for analyzing multivariate repeated measures designs. Multivariate Behavioral Research, 36(1):01–27.
 Vallejo Seco et al. (2007) Vallejo Seco, G., Gras, J. A., and Ato García, M. (2007). Comparative robustness of recent methods for analyzing multivariate repeated measures designs. Educational and Psychological Measurement, 67(3):410–432.
 Xiao and Zhang (2016) Xiao, S. and Zhang, J.T. (2016). Modified Tests for Heteroscedastic TwoWay MANOVA. Journal of Advanced Statistics, 1(1):1.
 Yanagihara and Yuan (2005) Yanagihara, H. and Yuan, K.H. (2005). Three approximate solutions to the multivariate Behrens–Fisher problem. Communications in Statistics  Simulation and Computation, 34(4):975–988.
 Yao (1965) Yao, Y. (1965). An approximate degrees of freedom solution to the multivariate Behrens Fisher problem. Biometrika, 52(1/2):139–147.
 Zhang (2011) Zhang, J.T. (2011). Twoway MANOVA with unequal cell sizes and unequal cell covariance matrices. Technometrics, 53(4):426–439.
Comments
There are no comments yet.