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?
 A: 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.
A: 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.
