Notes on Testing for Multivariate Normality
It is necessary to check the assumption of multivariate normality prior to applying procedures which have the assumption of MVN, and common variance. For example, prior to applying MANOVA, procedures such as PCA or linear regression.
There are several methods available for checking for multivariate normality, which include visual inspection of Mahalanobis distance values for each observation from the expected chisq quantile, as well as performing tests available such as the Mardia, Henze-Zirkler and Royston tests which each use different measures to test for multivariate normality.
Example Plotting Mahalanobis Distance over Chisq Quantiles
One visual method of inspecting multivariate data for normality is to plot the mahalanobis distances obtained over the quantiles of the chisq distribution with 2 degrees of freedom.
The library implements the Mahalanobis distance measure via:
\[ dist(X) = \left\{ (X-\bar{X})' S_{X,X}^{-1} (X-\bar{X}) \right\}^{1/2} \]
Where \(S_{X,X}^{-1}\) is the inverse of the estimated variance-covariance matrix and \(\bar{X}\) is the estimated column mean vector for the matrix X.
The visual assessment can be demonstrated on two data sets. The first example is a subset of the iris data set, subsetted on the species “virginica”, in order to generate an example of a multivariate normal data sample.
// imports
import java.io.File
import scala.collection.mutable
import au.id.cxd.math.data.CsvReader
import au.id.cxd.math.data.MatrixReader
import au.id.cxd.math.function.transform.StandardisedNormalisation
import au.id.cxd.math.probability.analysis._
import breeze.linalg.{DenseMatrix, eigSym, inv, svd}
import au.id.cxd.math.function.distance.MahalanobisDistance
import au.id.cxd.math.probability.continuous.ChiSquare
val fileName:String = s"$basePath/data/iris_virginica.csv"
val mat = MatrixReader.readFileAt(fileName)
// extract the columns of interest.
val data = mat(::, 0 to 3)Once the data is loaded we then obtain the distance measures for each observation.
// get the distance measures for each example.
val distance = MahalanobisDistance(data)
// sort distance from lowest to highest.
val distSorted = distance.toArray.sorted
// chisquare quantiles
val n = distSorted.length
val x = (for (i <- 1.0 to n.toDouble by 1.0) yield (i-0.5)/n.toDouble).toArray
val df = data.cols
val chisq = ChiSquare(df)
val quantiles = x.map(chisq.invcdf(_))Now we will produce a plot of the ordered distance measure over the chisq quantiles.
distances <- sapply(distances, function(d) d^2)
m <- lm(quantiles ~ distances)
intercept <- m$coefficients[1]
slope <- m$coefficients[2]
plot(distances, quantiles, col="blue", main="QQPlot of squared mahalanobis distance over Chisq Quantiles df=4")
abline(0,1, col="red")
The plot above suggests the data is close to multivariate normal, since the squared mahalanobis distance closely follows the chisquared distribution of the same degree of freedom, with 4 attributes.
We can also repeat this process using R for comparison.
data1 <- read.csv("../data/iris_virginica.csv", header=TRUE)
data <- (data1[,1:4])
mu <- colMeans(data)
sigma <- cov(data)
dist <- mahalanobis(data, mu, sigma)
df <- ncol(data)
n <- length(dist)
u <- ((1:n)-0.5)/n
p <- qchisq(u,df)
distsorted <- sort(dist)
plot(distsorted,p,
col="blue",
main="QQ Plot of mahalnobis distance v chisq quantiles")
abline(0,1, col="red")
A QQPlot for non-MVN distributed data
The mandible data set is not multivariate normally distributed, as the groups do not share common variance (as shown in the manova notes).
The qqplot can be displayed in a similar manner using the Mahalanobis distance measures.
val fileName2:String = s"$basePath/data/test_mandible_data.csv"
val dataRange:Seq[Int] = 2 to 10
val mat2 = MatrixReader.readFileAt(fileName2)
val m2 = mat2(::, dataRange).toDenseMatrix
val data2 = StandardisedNormalisation().transform(m2)
val distance2 = MahalanobisDistance(data2)
// sort distance from lowest to highest.
val distSorted2 = distance2.toArray.sorted
// chisquare quantiles
val n2 = distSorted2.length
val x2 = (for (i <- 1.0 to n2.toDouble by 1.0) yield (i-0.5)/n2.toDouble).toArray
val df2 = data2.cols
val chisq2 = ChiSquare(df2)
val quantiles2 = x2.map(chisq2.invcdf(_))Displaying the plot below
distances2 <- sapply(distances2, function(d) d^2)
plot(distances2, quantiles2, col="blue", main="QQPlot of mandible data")
abline(0,1, col="red")
It is apparent that there are a small number of examples whose distances deviate from the expected quantiles. It may be interesting to identify which observations they are and when the segments (species of dog in this case) are known before hand identify whether there may be some commonality of segment or not for those cases. It appears that the other samples may correspond closely to observations from an MVN distribution, it is possible that the four examples that deviate largely could be outliers or particularly interesting observations.
Tests for Multivariate Normality.
There are a number of procedures available for testing multivariate normality. The procedures are summarised in the article by Kormaz, Goksuluk and Zarasiz [1].
The procedures include
- The Mardia Test
- The Henze-Zirkler Test
- The Royston Test
When inspecting data for multivariate normality, it is advisable to use the combination of tests and visualisation, as well as using more than one test in order to determine whether the procedures agree.
The Mardia Test makes use of statistics for skew and kurtosis.
Like the Mahalanobis distance measure \(m\), the kurtosis and skew statistics are derived from the matrix
\[ m^2 = R = (X-\bar{X})'S^{-1}(X-\bar{X}) \]
Where we have the estimate of the parameters for the MVN distribution, the mean vector parameter \(\bar{X}\) and the estimate of the variance covariance matrix \(S\).
The measures of skewness \(b_{1,p}\) and kurtosis \(b_{2,p}\) are calculated as
\[ b_{1,p} = \frac{1}{n^2} \sum_{i,j=1}^n r^3_{ij} \] and \[ b_{2,p} = \frac{1}{n}\sum_{i=1}^n r^2_{ii} \] where \(n\) is the number of observations in the sample \(X\)
For large \(n\) the test statistic for skewness is derived as \[ z_1 = \frac{n}{6}b_{1,p} \] which has a \(\chi^2\) distribution with df \[ df = p(p+1)(p+2)/6 \]
and the statistic for kurtosis \(z_2\) is \[ z_2 = \sqrt{n} \frac{b_{2,p} - p(p+2) } { \sqrt{8 p(p+2) } } \] which is distributed as \(N(0,1)\). It is a standardisation of \(b_{2,p}\) with \(\mu=p(p+2)\) and \(\sigma = \sqrt{ 8p(p+2)/n}\)
In the case of small sample sizes less then 20 observations, a correction factor \(c\) is given for \(z_1\) \[ c = \frac{(n+1)(n+3)(k+1)}{n(n+1)(k+1)-6} \] then \(z1 = \frac{nc}{6}b_{1,p}\).
The null hypothesis asserts that the data is MVN distributed. Two tests are performed using the statistics and associated distributions. In order to reject the null hypothesis, either of the tests must disagree, both tests must agree in order to accept the MVN hypothesis.
import au.id.cxd.math.probability.analysis.MardiaTest
// test with standardised data
val X = StandardisedNormalisation().transform(data)
val test = MardiaTest(0.05, X)Printing the result:
println(test.toString)##
## Mardia Test at alpha = 0.05
## Skewness: 2.8444717653387688
## Skew Statistic: 23.703931377823075
## Skew DF: 20.0
## Skew P-Value: 0.2534434745725782
## Skew Critical Value: 31.37599076506435
## Skew Test Reject H0:false
##
## Kurtosis: 21.973352110762885
## Kurtosis Statistic: -1.0342194201915607
## Kurtosis P-Value: 0.8494832013464673
## Kurtosis Critical Value: 1.9599639523015855
## Kurtosis Test Reject H0: false
##
## Result: Cannot reject null hypothesis, sample is multivariate normal
## As both tests do not reject the null hypothesis, the outcome suggests that the set of observations are from a multivariate normal distribution.
The Mardia Test from the MVN package is demonstrated below for the iris subset of data from above. The results in the mardia test require at least one of the kurtosis and skew test results to disagree in order to for the null hypothesis to be rejected.
require(MVN)
mvn(data, mvnTest="mardia", multivariatePlot="qq")
## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 25.1850115362466 0.194444483140265 YES
## 2 Mardia Kurtosis -0.57186635893429 0.567412516528727 YES
## 3 MVN <NA> <NA> YES
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk Sepal.Length 0.9778 0.4647 YES
## 2 Shapiro-Wilk Sepal.Width 0.9741 0.3380 YES
## 3 Shapiro-Wilk Petal.Length 0.9660 0.1585 YES
## 4 Shapiro-Wilk Petal.Width 0.9476 0.0273 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew
## Sepal.Length 50 5.936 0.5161711 5.90 4.9 7.0 5.600 6.3 0.09913926
## Sepal.Width 50 2.770 0.3137983 2.80 2.0 3.4 2.525 3.0 -0.34136443
## Petal.Length 50 4.260 0.4699110 4.35 3.0 5.1 4.000 4.6 -0.57060243
## Petal.Width 50 1.326 0.1977527 1.30 1.0 1.8 1.200 1.5 -0.02933377
## Kurtosis
## Sepal.Length -0.6939138
## Sepal.Width -0.5493203
## Petal.Length -0.1902555
## Petal.Width -0.5873144An example of testing non-MVN data.
Similarly we can observe the result for the non-multivariate data for the mandible dataset.
// test with standardised data
val test2 = MardiaTest(0.05, data2)Printing the results;
println(test2)##
## Mardia Test at alpha = 0.05
## Skewness: 23.979566313845808
## Skew Statistic: 307.73776769435455
## Skew DF: 165.0
## Skew P-Value: 1.0708556263949731E-10
## Skew Critical Value: 195.9633891051611
## Skew Test Reject H0:true
##
## Kurtosis: 108.82316943260622
## Kurtosis Statistic: 3.062911205753588
## Kurtosis P-Value: 0.0010959754505143193
## Kurtosis Critical Value: 1.9599639523015855
## Kurtosis Test Reject H0: true
##
## Result: Reject null hypothesis, sample is not multivariate normal
## We can see that the test rejects the null hypothesis of MVN for this particulare dataset.
For comparison in R as below,
data2 <- read.csv("../data/test_mandible_data.csv", header=TRUE)
mvn(data2[,1:9], mvnTest="mardia")## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 276.023492238367 1.32921080879424e-07 NO
## 2 Mardia Kurtosis 1.55269100090827 0.12049697326799 YES
## 3 MVN <NA> <NA> NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk Case 0.9628 0.0237 NO
## 2 Shapiro-Wilk Group 0.8978 <0.001 NO
## 3 Shapiro-Wilk X1 0.9208 1e-04 NO
## 4 Shapiro-Wilk X2 0.9824 0.3672 YES
## 5 Shapiro-Wilk X3 0.9359 8e-04 NO
## 6 Shapiro-Wilk X4 0.9538 0.007 NO
## 7 Shapiro-Wilk X5 0.9116 1e-04 NO
## 8 Shapiro-Wilk X6 0.9861 0.5705 YES
## 9 Shapiro-Wilk X7 0.8588 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew
## Case 77 8.558442 5.014436 8.0 1.0 20.0 4.0 12.0 0.28017061
## Group 77 2.766234 1.326810 3.0 1.0 5.0 2.0 4.0 0.22691891
## X1 77 128.974026 17.501860 125.0 105.0 177.0 114.0 137.0 0.82425183
## X2 77 9.961039 1.403019 10.0 7.2 13.4 8.7 10.9 0.03445269
## X3 77 21.948052 3.583207 22.0 17.0 32.0 19.0 25.0 0.67695177
## X4 77 21.493506 3.378038 22.0 15.0 28.0 18.0 24.0 -0.14798176
## X5 77 20.493506 2.490103 20.0 17.0 27.0 19.0 22.0 0.75957626
## X6 77 8.000000 1.023667 7.9 6.0 10.5 7.1 8.7 0.14506482
## X7 77 32.519481 4.172625 31.0 26.0 43.0 30.0 33.0 1.08564435
## Kurtosis
## Case -0.92889918
## Group -1.13911095
## X1 -0.06450772
## X2 -0.75431906
## X3 -0.23710896
## X4 -1.13282761
## X5 -0.39563109
## X6 -0.64113198
## X7 0.23084116reinforces the conclusion that the dataset is not multivariate normal.
References
[1] Korkmaz, S, Goksuluk, D and Zararsiz, G (2014 ). MVN: An r package for assessing multivariate normality. The R Journal. 6 151–62. http://journal.r-project.org/archive/2014-2/korkmaz-goksuluk-zararsiz.pdf