Suppose in R, I have a data set


and I want to check for multivariate normality of the first four columns ("Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width") using the rQ statistic. How can one do this?

I understand how to check if each of the marginals are normal (for example, ppcc function in R will do this), but I dont see how I can use rQ to check for multivariate normality. Supposing one of the ppcc calls fails marginal normality, then we can reject multivariate normality; but even if all four of the marginals are normal (checked by rQ) this does not imply multivariate normality.

How do I proceed?

  • $\begingroup$ Could you please define what is the rQ statistic? $\endgroup$ – Lucas Farias Feb 12 at 4:19

I think it's also important to flesh out the underlying statistical motivation for the various multivariate tests already mentioned - beyond just how to carry it out in a statistical software.

Let's assume you have a multivariate Normal distribution

$$\boldsymbol{X}\sim N_{d}(\boldsymbol{\mu},\boldsymbol{\Sigma})$$



So we have that the above is Chi-squared distributed with $d$ degrees of freedom. We can then exploit this fact and test for joint-normality. We can calculate


for $i=1,\ldots,n$, where $\hat{\boldsymbol{X}}$ and $S$ are your estimates for $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, respectively. Technically speaking, because of the use of the above estimates, the $D_{i}$ are not independent of each other even if the original $X_{i}$ are independent. $D_{i}$ is generally known as the Mahalanobis distance. Thus, under the null hypothesis of multivariate Normality, we have that


However, for large $n$, we have that it is well approximated by a Chi-squared distribution as above. Given this, we can construct QQ-plots of the $D_{i}^{2}$ against that of a Chi-squared (or Beta) distribution. An example is shown below using R:

dd=dif %*% solve(S) %*% t(dif)
plot(d,chi2q,pch=20,main="",xlab="Mahalanobis distance",ylab="Chi-squared quantile",col="blue")

Mahalanobis distance

Mardia's test of multivariate normality is based on multivariate measures of skewness and kurtosis and compares these with theoretical reference distributions (much like how the Jarque-Bera test assesses univariate normality).

If we define the following measures

$$b_{d}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}D_{i,j}^{3}$$ $$k_{d}=\frac{1}{n}\sum_{i=1}^{n}D_{i}^{4}$$

Under the null hypothesis of multivariate normality we have that as $n\rightarrow\infty$

$$\frac{1}{6}nb_{d}\sim\chi^{2}_{d(d+1)(d+2)/6}$$ $$\frac{k_{d}-d(d+2)}{\sqrt{8d(d+2)/n}}\sim N(0,1)$$

Generally, tests of the skewness and kurtosis are performed separately but there are joint tests too.

  • $\begingroup$ Nice explanation! $\endgroup$ – Lucas Farias Feb 12 at 14:15

I believe the package MVN is what you're looking for. Following from your example dataset, which is actually the one they use in their vignette:

# load MVN package

# load Iris data

# setosa subset of the Iris data
setosa <- iris[1:50, 1:4]

result <- mvn(data = setosa, mvnTest = "mardia")

             Test        Statistic           p value Result
1 Mardia Skewness 25.6643445196298 0.177185884467652    YES
2 Mardia Kurtosis 1.29499223711605 0.195322907441935    YES
3             MVN             <NA>              <NA>    YES

There is a variety of tests implemented: Mardia’s test, Henze-Zirkler’s test, Royston’s test, Doornik-Hansen’s test and the E-statistic test.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.