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

01/24/2018 ∙ by Sarah Friedrich, et al. ∙ Universität Ulm 0

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.



There are no comments yet.


page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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:

  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.

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

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 one- and 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, 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 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); 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); Livacic-Rojas et al. (2017); Friedrich and Pauly (2017)

and the references cited therein. Here, we focus on statistical methods that are valid in the multivariate Behrens-Fisher situation—equal covariance matrices across the groups is not assumed—and provide accurate inferential results in terms of p-value estimates and confidence intervals for the parameters of interest. 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 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 Behrens-Fisher 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

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 user-friendly application not only provides the -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 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 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 -distribution for -value computations, and

  • the modified ANOVA-type statistic (MATS) proposed by Friedrich and Pauly (2017) using a parametric or wild bootstrap procedure for -value computations.

In addition to multivariate group-wise effects, the RM function also allows to test hypotheses formulated across sub-plot factors. The implemented test statistics are

  • the ANOVA-type statistic (ATS) using an -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 -distribution (Brunner, 2001), the permutation technique proposed in Friedrich et al. (2017a) as well as a parametric (Bathke et al., 2016) and a wild bootstrap approach for -value estimation.

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


Here, denotes the experimental unit or subject in group . Note, that a higher-way factorial structure on the groups/whole-plots or sub-plots can be achieved by sub-indexing the indices (group/whole-plot) or (sub-plot) 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 WTS-type 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 set-up Konietschke et al. (2015) propose

a Wald-type statistic (WTS)


for testing , where , , and denotes the Moore-Penrose inverse of the matrix . Since its asymptotic null distribution provides a poor finite sample approximation, they propose the following asymptotic model-based 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 type-1 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 well-established ANOVA-type statistic (ATS)


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 Moore-Penrose 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)


Here, the involved diagonal matrix of the empirical variances of component in group , deduces an invariance under component-wise 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 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 in model (1) in the following ways:

  • One-Way (): Writing we have with

    and obtain the null hypothesis of ’no group’ or ’factor

    ’ effect as

    In case of 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 (): 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 Two-Way (): A fixed subcategory within factor can be introduced via the model with Here, the hypotheses of no main effect or no sub-category main effect can be written as

    with and .

    We only implemented balanced designs, i.e., for all . 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 , where even hypotheses about sub-plots can be formulated. We exemplify this for profile analyses in the special case of a one-sample RM design with

    as well as for a two-sample 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 sub-plot 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 higher-way 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 = ”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 sub-plot 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 whole-plot and two sub-plot 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 sub-plot factors while the treatment (Verum vs. Placebo) is a whole-plot 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

Wald-Type Statistic (WTS): Test statistic df p-value 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

ANOVA-Type Statistic (ATS): Test statistic df1 df2 p-value 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

p-values 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 Wald-type 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 sub-plot and two whole-plot 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 whole-plot factors considered were sex (men vs. women), diagnosis (AD vs. MCI vs. SCC), and age ( vs. years). Additionally, the sub-plot 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

Wald-Type Statistic (WTS): Test statistic df p-value 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

ANOVA-Type Statistic (ATS): Test statistic df1 df2 p-value 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

p-values 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 whole-plot factors sex and diagnosis, while none of the sub-plot 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 higher-way 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. 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)

The resulting plots are displayed in Figure 1 and Figure 2, respectively.

Figure 1: Plot for factor ”sex” in the RM model of the EEG data example.
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.

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 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 (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 sub-plot 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

Wald-Type Statistic (WTS): Test statistic df p-value sex 12.604 6 0.050 diagnosis 55.158 12 0.000 sex:diagnosis 9.790 12 0.634

modified ANOVA-Type Statistic (MATS): Test statistic sex 45.263 diagnosis 194.165 sex:diagnosis 18.401

p-values 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 resampling-based -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

Wald-Type Statistic (WTS): Test statistic df p-value sex 12.604 6 0.050 diagnosis 55.158 12 0.000 sex:diagnosis 9.790 12 0.634

modified ANOVA-Type Statistic (MATS): Test statistic sex 45.263 diagnosis 194.165 sex:diagnosis 18.401

p-values 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

Wald-Type Statistic (WTS): Test statistic df p-value 51.584 2.000 0.000

modified ANOVA-Type Statistic (MATS): [,1] [1,] 69.882

p-values 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 two-dimensional outcomes as in this example, the confidence ellipsoid can also be plotted, see Figure


Center: [,1] [1,] 256.792 [2,] -39.369

Scale: [1] 10.852716 2.736354

Eigenvectors: [,1] [,2] [1,] -1 0 [2,] 0 -1

Figure 3: Plot of the confidence region for factor location in the water example from package HSAUR.

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

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

Wald-Type Statistic (WTS): Test statistic df p-value season 6.999 2 0.030 season:site 16.621 8 0.034

modified ANOVA-Type Statistic (MATS): Test statistic season 12.296 season:site 15.064

p-values 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

Wald-Type Statistic (WTS): Test statistic df p-value season 6.999 2 0.030 season:site 16.621 8 0.034

modified ANOVA-Type Statistic (MATS): Test statistic season 12.296 season:site 15.064

p-values 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 sub-plot 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.

R¿ library(”MANOVA.RM”) R¿ GUI.RM()

Figure 4: The GUI for tests in repeated measures designs (upper panel) and multivariate data (lower panel): The user can specify the data location and the formula as well as the resampling approach.
Figure 5: Graphical user interfaces for plotting: 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.

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 ANOVA-type 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 sub-plot effects in RM designs with possibly complex factorial structures on both, whole- and sub-plots.

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.


The work of Sarah Friedrich and Markus Pauly was supported by the German Research Foundation project DFG-PA 2409/3-1.


  • Anderson (2001) Anderson, M. J. (2001). A new method for non-parametric 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). Box-type 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 Behrens-Fisher 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.3-7.
  • 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 non-normality. Test, 17(1):120–137.
  • Harrar and Bathke (2012) Harrar, S. W. and Bathke, A. C. (2012). A modified two-factor 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. R-News, 6(4):8–12.
  • Livacic-Rojas et al. (2010) Livacic-Rojas, 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.
  • Livacic-Rojas et al. (2017) Livacic-Rojas, P., Vallejo, G., Fernández, P., and Tuero-Herrero, E. (2017). Power of modified brown-forsythe and mixed-model approaches in split-plot 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 group-level 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 Behrens-Fisher problem. Communications in Statistics-Theory 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 Two-Way 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). Two-way MANOVA with unequal cell sizes and unequal cell covariance matrices. Technometrics, 53(4):426–439.