# Using $r_Q$ statistic to check multivariate normality

Suppose in R, I have a data set

data(iris)


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?

• Could you please define what is the rQ statistic? – 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})$$

then

$$(\boldsymbol{X}-\boldsymbol{\mu})\boldsymbol{\Sigma}^{-1}(\boldsymbol{X}-\boldsymbol{\mu})\sim\chi^{2}_{d}$$

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

$$D_{i}^{2}=(\boldsymbol{X}_{i}-\hat{\boldsymbol{X}})'S^{-1}(\boldsymbol{X}_{i}-\hat{\boldsymbol{X}})$$

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

$$n(n-1)^{-2}D_{i}^{2}\sim\text{Beta}(d/2,(n-d-1)/2)$$

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:

n=dim(data)[1]
p=dim(data)[2]
S=cov(data)
dif=scale(data,scale=FALSE)
dd=dif %*% solve(S) %*% t(dif)
d=diag(dd)
r=rank(d)
chi2q=qchisq((r-0.5)/n,p)
plot(d,chi2q,pch=20,main="",xlab="Mahalanobis distance",ylab="Chi-squared quantile",col="blue")
abline(0,1,lwd=2,col="red")


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.

• Nice explanation! – 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
library(MVN)

# load Iris data
data(iris)

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

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

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.