MVN : An R Package for Assessing Multivariate Normality

Assessing the assumption of multivariate normality is required by many parametric multivariate statistical methods, such as MANOVA, linear discriminant analysis, principal component analysis, canonical correlation, etc. It is important to assess multivariate normality in order to proceed with such statistical methods. There are many analytical methods proposed for checking multivariate normality. However, deciding which method to use is a challenging process, since each method may give different results under certain conditions. Hence, we may say that there is no best method, which is valid under any condition, for normality checking. In addition to numerical results, it is very useful to use graphical methods to decide on multivariate normality. Combining the numerical results from several methods with graphical approaches can be useful and provide more reliable decisions. Here, we present an R package, MVN, to assess multivariate normality. It contains the three most widely used multivariate normality tests, including Mardia’s, Henze-Zirkler’s and Royston’s, and graphical approaches, including chi-square Q-Q, perspective and contour plots. It also includes two multivariate outlier detection methods, which are based on robust Mahalanobis distances. Moreover, this package offers functions to check the univariate normality of marginal distributions through both tests and plots. Furthermore, especially for non-R users, we provide a user-friendly web application of the package. This application is available at http://www.biosoft.hacettepe.edu.tr/MVN/.


Introduction
Many multivariate statistical analysis methods, such as MANOVA and linear discriminant analysis (MASS, [1]), principal component analysis (FactoMineR, [2], psych, [3]), canonical correlation (CCA, [4]), etc., require multivariate normality (MVN) assumption.If the data are multivariate normal (exactly or approximately), such multivariate methods provide more reliable results.The performance of these methods dramatically decreases if the data are not multivariate normal.Hence, researchers should check whether data are multivariate normal or not before continuing with such parametric multivariate analyses.
Many statistical tests and graphical approaches are available to check the multivariate normality assumption.Burdenski (2000) reviewed several statistical and practical approaches, including the Q-Q plot, box-plot, stem and leaf plot, Shapiro-Wilk and Kolmogorov-Smirnov tests to evaluate the univariate normality, contour and perspective plots for assessing bivariate normality, and the chisquare Q-Q plot to check the multivariate normality [5].The author demonstrated each procedure using the real data from [6].Ramzan et al. (2013) reviewed numerous graphical methods for assessing both univariate and multivariate normality and showed their use in a real life problem to check the MVN using chi-square and beta Q-Q plots [7].Holgersson (2006) stated the importance of graphical procedures and presented a simple graphical tool, which is based on the scatter plot of two correlated variables to assess whether the data belong to a multivariate normal distribution or not [8].Svantesson and Wallace (2003) applied Royston's and Henze-Zirkler's tests to multipleinput multiple-output data to test MVN [9].According to the review by Mecklin and Mundfrom (2005), more than fifty statistical methods are available for testing MVN [10].They conducted a comprehensive simulation study based on type I and type II error and concluded that no single test excelled in all situations.The authors suggested using Henze-Zirkler's and Royston's tests among others for assessing MVN because of their good type I error control and power.Moreover, to diagnose the reason for deviation from multivariate normality, the authors suggested the use of Mardia's multivariate skewness and kurtosis statistics test as well as graphical approaches such as the chi-square Q-Q plot.Deciding which test to use can be a daunting task for researchers (mostly for non-statisticians) and it is very useful to perform several tests and examine the graphical methods simultaneously.Although there are a number of studies describing multifarious approaches, there is no single easy-to-use, up-to-date and comprehensive tool to apply various statistical tests and graphical methods together at present.
In this vignette, we introduce an R package, MVN, which implements the three most widely used MVN tests, including Mardia's, Henze-Zirkler's, and Royston's [11].In addition to statistical tests, the MVN also provides some graphical approaches such as chi-square Q-Q, perspective and contour plots.Moreover, this package includes two multivariate outlier detection methods, which are based on Mahalanobis distance.In addition to multivariate normality, users can also check univariate normality tests and plots to diagnose the deviation from normality via package version 3.7 and later.Firstly, we discuss the theoretical background on the corresponding MVN tests.Secondly, two illustrative examples are presented in order to demonstrate the applicability of the package.Finally, we present a newly developed web interface of the MVN, which can be especially handy for non-R users.The R version of the MVN is publicly available in the Comprehensive R Archive Network (CRAN, http://CRAN.R-project.org/package=MVN).

Henze-Zirkler's MVN test
The Henze-Zirkler's test is based on a non-negative functional distance that measures the distance between two distribution functions.If data are distributed as multivariate normal, the test statistic is approximately log-normally distributed.First, the mean, variance and smoothness parameter are calculated.Then, the mean and the variance are log-normalized and the p-value is estimated [14][15][16][17][18]. The test statistic of Henze-Zirkler's multivariate normality test is given in equation 2.
where p : number of variables From equation 2, D i gives the squared Mahalanobis distance of i th observation to the centroid and D ij gives the Mahalanobis distance between i th and j th observations.If data are multivariate normal, the test statistic (HZ) is approximately log-normally distributed with mean µ and variance σ 2 as given below: where a = 1 + 2β 2 and w β = (1 + β 2 )(1 + 3β 2 ).Hence, the log-normalized mean and variance of the HZ statistic can be defined as follows: By using the log-normal distribution parameters, µ and σ, we can test the significance of multivariate normality.The Wald test statistic for multivariate normality is given in equation 4.

Royston's MVN test
Royston's test uses the Shapiro-Wilk/Shapiro-Francia statistic to test multivariate normality.If kurtosis of the data is greater than 3, then it uses the Shapiro-Francia test for leptokurtic distributions, otherwise it uses the Shapiro-Wilk test for platykurtic distributions [10,15,[19][20][21][22][23].Let W j be the Shapiro-Wilk/Shapiro-Francia test statistic for the j th variable ( j = 1, 2, . . ., p) and Z j be the values obtained from the normality transformation proposed by [22].
x = log(n) and As seen from equation 5, x and w j 's change with the sample size (n).By using equation 5, transformed values of each random variable can be obtained from equation 6.
where γ, µ and σ are derived from the polynomial approximations given in equation 7. The polynomial coefficients are provided by [22] for different sample sizes.
The Royston's test statistic for multivariate normality as follows: where e is the equivalent degrees of freedom (edf) and Φ(.) is the cumulative distribution function for standard normal distribution such that, As seen from equation 9, another extra term c has to be calculated in order to continue with the statistical significance of Royston's test statistic given in equation 8. Let R be the correlation matrix and r ij be the correlation between i th and j th variables.Then, the extra term c can be found by using equation 10. where with the boundaries of g(.) as g(0, n) = 0 and g(1, n) = 1.The function g(.) is defined as follows: The unknown parameters, µ, λ and ν were estimated from a simulation study conducted by [24].He found µ = 0.715 and λ = 5 for sample size 10 ≤ n ≤ 2000 and ν is a cubic function which can be obtained as follows: where x = log(n).

Implementation of MVN package
The MVN package contains several functions in the S4 class.The data to be analyzed should be given in the "data.frame"or "matrix" class.In this example, we will work with the famous Iris data set.These data are from a multivariate data set introduced by Fisher (1936) as an application of linear discriminant analysis [25].It is also called Anderson's Iris data set because Edgar Anderson collected the data to measure the morphologic variation of Iris flowers of three related species [26].First of all, the MVN library should be loaded in order to use related functions.

# load MVN package library(MVN)
Similarly, Iris data can be loaded from the R database by using the following R code: The Iris data set consists of 150 samples from each of the three species of Iris including setosa, virginica and versicolor.For each sample, four variables were measured including the length and width of the sepals and petals, in centimeters.

Chi-square Q-Q plot
One can clearly see that different MVN tests may come up with different results.MVN assumption was rejected by Henze-Zirkler's and Royston's tests; however, it was not rejected by Mardia's test at a significance level of 0.05.In such cases, examining MVN plots along with hypothesis tests can be quite useful in order to reach a more reliable decision.
The Q-Q plot, where "Q" stands for quantile, is a widely used graphical approach to evaluate the agreement between two probability distributions.Each axis refers to the quantiles of probability distributions to be compared, where one of the axes indicates theoretical quantiles (hypothesized quantiles) and the other indicates the observed quantiles.If the observed data fit hypothesized distribution, the points in the Q-Q plot will approximately lie on the line y = x.
MVN has the ability to create three multivariate plots.One may use the qqplot = TRUE option in the mardiaTest, hzTest and roystonTest functions to create a chi-square Q-Q plot.We can create this plot for the setosa data set to see whether there are any deviations from multivariate normality.Figure 1 shows the chi-square Q-Q plot of the first 50 rows of Iris data, which are setosa flowers.It can be seen from Figure 1 that there are some deviations from the straight line and this indicates possible departures from a multivariate normal distribution.As a result, we can conclude that this data set does not satisfy MVN assumption based on the fact that the two test results are against it and the chi-square Q-Q plot indicates departures from multivariate normal distribution.

Univariate plots and tests
As noted by several authors [5,27,28], if data have a multivariate normal distribution, then, each of the variables has a univariate normal distribution; but the opposite does not have to be true.Hence, checking univariate plots and tests could be very useful to diagnose the reason for deviation from MVN.We can check this assumption through uniPlot and uniNorm functions from the package.The uniPlot function is used to create univariate plots, such as Q-Q plots (Figure 2a), histograms with normal curves (Figure 2b), box-plots and scatterplot matrices.
uniPlot(setosa, type = "qqplot") # creates univariate Q-Q plots uniPlot(setosa, type = "histogram") # creates univariate histograms  As seen from Figure 2, Petal.Width has a right-skewed distribution whereas other variables have approximately normal distributions.Thus, we can conclude that problems with multivariate normality arise from the skewed distribution of Petal.Width.In addition to the univariate plots, one can also perform univariate normality tests using the uniNorm function.It provides several widely used univariate normality tests, including Shapiro-Wilk, Cramer-von Mises, Lilliefors and Anderson-Darling.For example, the following code chunk is used to perform the Shapiro-Wilk's normality test on each variable and it also displays descriptive statistics including mean, standard deviation, median, minimum, maximum, 25th and 75th percentiles, skewness and kurtosis: uniNorm(setosa, type = "SW", desc = TRUE) Example II: Whilst the Q-Q plot is a general approach for assessing MVN in all types of numerical multivariate datasets, perspective and contour plots can only be used for bivariate data.To demonstrate the applicability of these two approaches, we will use a subset of Iris data, named setosa2, including the sepal length and sepal width variables of the setosa species.

Perspective and contour plots
Univariate normal marginal densities are a necessary but not a sufficient condition for MVN.Hence, in addition to univariate plots, creating perspective and contour plots will be useful.The perspective plot is an extension of the univariate probability distribution curve into a 3•dimensional probability distribution surface related with bivariate distributions.It also gives information about where data are gathered and how two variables are correlated with each other.It consists of three dimensions where two dimensions refer to the values of the two variables and the third dimension, which is likely in univariate cases, is the value of the multivariate normal probability density function.Another alternative graph, which is called the "contour plot", involves the projection of the perspective plot into a 2•dimensional space and this can be used for checking multivariate normality assumption.For bivariate normally distributed data, we expect to obtain a three-dimensional bell-shaped graph from the perspective plot.Similarly, in the contour plot, we can observe a similar pattern.
To construct a perspective and contour plot for Example 2, we can use the mvnPlot function in the MVN.This function requires an object in the "MVN" class that is one of the results from MVN functions.In the following codes, the object from hzTest is used for the perspective plot given in Figure 3a.It is also possible to create a contour plot of the data.Contour graphs are very useful since they give information about normality and correlation at the same time.Figure 3b shows the contour plot of setosa flowers.As can be seen from the graph, this is simply a top view of the perspective plot where the third dimension is represented with ellipsoid contour lines.From this graph, we can say that there is a positive correlation among the sepal measures of flowers since the contour lines lie around the main diagonal.If the correlation were zero, the contour lines would be circular rather than ellipsoid.Since neither the univariate plots in Figure 2 nor the multivariate plots in Figure 3 show any significant deviation from MVN, we can now perform the MVN tests to evaluate the statistical significance of bivariate normal distribution of the setosa2 data set.Figures 3a and 3b were drawn using a pre-defined graphical option by the authors.However, users may change these options by setting function entry to default = FALSE.If the default is FALSE, optional arguments from the plot, persp and contour functions may be introduced to the corresponding graphs.

Multivariate outliers
Multivariate outliers are the common reason for violating MVN assumption.In other words, MVN assumption requires the absence of multivariate outliers.Thus, it is crucial to check whether the data have multivariate outliers, before starting to multivariate analysis.The MVN includes two multivariate outlier detection methods which are based on robust Mahalanobis distances (rMD(x)).Mahalanobis distance is a metric which calculates how far each observation is to the center of joint distribution, which can be thought of as the centroid in multivariate space.Robust distances are estimated from minimum covariance determinant estimators rather than the sample covariance [29].These two approaches, defined as Mahalanobis distance and adjusted Mahalanobis distance in the package, detect multivariate outliers as given below, Mahalanobis Distance: The mvOutlier function is used to detect multivariate outliers as given below.It also returns a new data set in which declared outliers are removed.Moreover, Q-Q plots can be created by setting qqplot = TRUE within mvOutlier for visual inspection of the possible outliers.For this example, we will use another subset of the Iris data, which is versicolor flowers, with the first three variables.

Web interface for the MVN package
The purpose of the package is to provide MVN tests along with graphical approaches for assessing MVN.Moreover, this package offers univariate tests and plots, and multivariate outlier detection for checking MVN assumptions through R.However, using R codes might be challenging for new R     users.Therefore, we also developed a user-friendly web application of MVN by using shiny1 [31].This web-tool, which is an interactive application, has all the features that the MVN package has.

Summary and further researches
As stated earlier, MVN is among the most crucial assumptions for most parametric multivariate statistical procedures.The power of these procedures is negatively affected if this assumption is not satisfied.Thus, before using any of the parametric multivariate statistical methods, MVN assumption should be tested first of all.Although there are many MVN tests, there is not a standard test for assessing this assumption.In our experience, researchers may choose Royston's test for data with a small sample size (n < 50) and Henze-Zirkler's test for a large sample size (n > 100).However, a more comprehensive simulation study is needed to provide more reliable inference.Instead of using just one test, it is suggested that using several tests simultaneously and examining some graphical representation of the data may be more appropriate.Currently, as we know, there is no such extensive tool to apply different statistical tests and graphical methods together.
In this vignette, we present the MVN package for multivariate normality checking.This package offers comprehensive flexibility for assessing MVN assumption.It contains the three most widely used MVN tests, including Mardia's, Henze-Zirkler's and Royston's.Moreover, researchers can create three MVN plots using this package, including the chi-square Q-Q plot for any data set and perspective and contour plots for bivariate data sets.Furthermore, since MVN requires univariate normality of each variable, users can check univariate normality assumption by using both univariate normality tests and plots with proper functions in the package.In the first example, different results on multivariate normality were achieved from the same data.When sepal and petal measures, i.e four variables, were considered, Mardia's test resulted in multivariate normality as well as Henze-Zirkler's test at the edge of type I error.However, Royston's test strongly rejected the null hypothesis in favor of non-normality.At this point, the only possible graphical approach is to use the chisquare Q-Q plot since there are more than two variables.The next step was to identify the cause of deviation from MVN by using univariate normality tests and plots.In the second example, all tests suggested bivariate normality, as did the graphical approaches.Although some tests can not reject null hypothesis, other tests may reject it.Hence, as stated earlier, selecting the appropriate MVN test dramatically changes the results and the final decision is ultimately the researcher's.
Currently, MVN works with several statistical tests and graphical approaches.It will continue to add new statistical approaches as they are developed.The package and the web-tool will be regularly updated based on these changes.

Figure 1 :
Figure 1: Chi-Square Q-Q plot for setosa data set.
Histograms with normal curves.

Figure 3 :
Figure 3: Perspective and contour plot for bivariate setosa2 data set.

:
-Mardia's estimation of multivariate skewness, i.e γ1,p given in equation 1, chi.skew: test statistic for multivariate skewness, p.value.skew: significance value of skewness statistic, g2p: Mardia's estimation of multivariate kurtosis, i.e γ2,p given in equation 1, z.kurtosis: test statistic for multivariate kurtosis, p.value.kurt: significance value of kurtosis statistic, chi.small.skew:test statistic for multivariate skewness with small sample correction, p.value.small: significance value of small sample skewness statistic.As seen from the results given above, both the skewness (γ 1,p = 3.0797, p = 0.1772) and kurtosis (γ 2,p = 26.5377,p = 0.1953) estimates indicate multivariate normality.Therefore, according to Mardia's MVN test, this data set follows a multivariate normal distribution.One may use the hzTest function in the MVN to perform the Henze-Zirkler's test.

Table 1 :
From the above results, we can see that all variables, except Petal.Width in the setosa data set, have univariate normal distributions at significance level 0.05.We can now drop Petal.With from setosa data and recheck the multivariate normality.MVN results are given in Table1.MVN test results (setosa without Petal.Width).to the three MVN test results in Table1, setosa without Petal.Width has a multivariate normal distribution at significance level 0.05.

Table 2 :
MVN test results (setosa with sepal measures).All three tests inTable 2 indicate that the data set satisfies bivariate normality assumption at the significance level 0.05.Moreover, the perspective and contour plots are in agreement with the test results and indicate approximate bivariate normality.