R code for multivariate normality testing
The R code below performs Mardia's test of multivariate kurtosis:
mardiatest <- function(x,alpha=0.05) {
# MARDIATEST Mardia multivariate normality test using skewness & kurtosis.
# MARDIATEST(X,ALPHA) performs Mardia's test of multivariate
# normality to determine if the null hypothesis of multivariate
# normality is a reasonable assumption regarding the population
# distributions of a random samples contained within the columns of X. X
# must be a N-values by M-samples array. The desired significance level,
# ALPHA, is an optional scalar input (default = 0.05).
#
# The function performs three tests: of the multivariate skewness; the
# multivariate skewness corrected for small samples; and the multivariate
# kurtosis. H is a 1-by-3 array indicating the results of the hypothesis
# tests:
# H(i) = 0 => Do not reject the null hypothesis at significance level
# ALPHA.
# H(i) = 1 => Reject the null hypothesis at significance level ALPHA.
#
# MARDIATEST(...) also returns the following fields:
# stats.Hs - logical scalar indicating whether to reject the null
# hypothesis that the multivariate skewness is consistent
# with a multivariate normal distribution.
# stats.Ps - the P value for the multivariate skewness.
# stats.Ms - the Mardia test statistic for the multivariate skewness.
# stats.CVs - the critical value for the multivariate skewness.
# stats.Hsc - logical scalar indicating whether to reject the null
# hypothesis that the multivariate skewness (corrected for
# small samples) is consistent with a multivariate normal
# distribution.
# stats.Psc - the P value for the multivariate skewness (corrected for
# small samples).
# stats.Msc - the Mardia test statistic for the multivariate skewness
# (corrected for small samples).
# stats.Hk - logical scalar indicating whether to reject the null
# hypothesis that the multivariate kurtosis is consistent
# with a multivariate normal distribution.
# stats.Pk - the P value for the multivariate kurtosis.
# stats.Mk - the Mardia test statistic for the multivariate kurtosis.
# stats.CVk - the critical value for the multivariate kurtosis.
#
# Notes:
# For multivariate data, tests of normality of the individual variables
# are not sufficent to determine multivariate normality. Even if every
# variable in a set is normally distributed, it is still possible that
# the combined distribution is not multivariate normal. Mardia's test is
# one means testing for multivariate normality. The test is based on
# independent tests of multivariate skewness and kurtosis. Data can be
# assumed to conform to multivariate normality only if the null
# hypothesis is not rejected for both tests.
#
# Example:
# >> X = [2.4 2.1 2.4; ...
# 3.5 1.8 3.9; ...
# 6.7 3.6 5.9; ...
# 5.3 3.3 6.1; ...
# 5.2 4.1 6.4; ...
# 3.2 2.7 4.0; ...
# 4.5 4.9 5.7; ...
# 3.9 4.7 4.7; ...
# 4.0 3.6 2.9; ...
# 5.7 5.5 6.2; ...
# 2.4 2.9 3.2; ...
# 2.7 2.6 4.1];
#
# >> [H stats] = mardiatest(X, 0.05)
#
# H =
# 0 0 1
#
# stats =
# Hs: 0
# Ps: 0.9660
# Ms: 3.5319
# CVs: 18.3070
# Hsc: 0
# Psc: 0.8918
# Msc: 4.9908
# Hk: 1
# Pk: 0.0452
# Mk: -1.6936
# CVk: 1.6449
#
# The magnitude of the Mardia kurtosis test statistic (Mk) is greater
# than critical value (CVk) for a 5% level test, so the hypothesis of
# multivariate normality is rejected.
#
#
# This program is based upon the MATCLAB version which is:
# Copyright 2006 David Graham, Loughborough University
# (D.J.Graham@lboro.ac.uk)
#
# This function is based on MSKEKUR by A. Trujillo-Ortiz and R.
# Hernandez-Walls. Modifications include:
# - additional checking of inputs
# - modification of help text for consistency with Matlab conventions
# - the output of the test results to variables
# - removal of the code sending test results as a string to the command
# window
# - removed call to 'eval'.
#
# Statements from the original code:
#
# Created by A. Trujillo-Ortiz and R. Hernandez-Walls
# Facultad de Ciencias Marinas
# Universidad Autonoma de Baja California
# Apdo. Postal 453
# Ensenada, Baja California
# Mexico
# atrujo@uabc.mx
#
# May 22, 2003.
#
# To cite this file, this would be an appropriate format:
# Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). Mskekur: Mardia's multivariate skewness
# and kurtosis coefficients and its hypotheses testing. A MATLAB file. [WWW document]. URL
# http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3519&objectType=FILE
#
# References:
#
# Mardia, K. V. (1970), Measures of multivariate skewnees and kurtosis with
# applications. Biometrika, 57(3):519-530.
# Mardia, K. V. (1974), Applications of some measures of multivariate skewness
# and kurtosis for testing normality and robustness studies. Sankhyâ A,
# 36:115-128
# Stevens, J. (1992), Applied Multivariate Statistics for Social Sciences. 2nd. ed.
# New-Jersey:Lawrance Erlbaum Associates Publishers. pp. 247-248.
# Check for validity of ALPHA
if (is.numeric(alpha) == FALSE || alpha>0.5 || alpha <0)
cat('Input ALPHA must be a scalar between 0 and 0.5.')
end
n <- nrow(x)
p <- ncol(x)
# Check for validity of X
if (p < 2)
cat('Input X must be an array with at least 2 columns.')
end
covxinv <- solve(cov(x))
dift <- matrix(0,nrow(x),ncol(x))
for (j in 1:ncol(x)) {
dift[,j] <- (x[,j] - mean(x[,j]))
}
D <- dift %*% covxinv %*% t(dift) # Mahalanobis Distances
b1p <- (sum(sum(D^3))) / n^2 # Multivariate skewness coefficient
b2p <- sum(diag(D^2)) / n # Multivariate kurtosis coefficient
k <- ((p+1)*(n+1)*(n+3)) / (n*(((n+1)*(p+1))-6)) # Small sample correction
v <- (p*(p+1)*(p+2)) / 6 # Degrees of freedom
g1c <- (n*b1p*k) / 6 # Skewness test statistic corrected for small sample (approximates to a chi-square distribution)
g1 <- (n*b1p) / 6 # Skewness test statistic (approximates to a chi-square distribution)
P1 <- 1 - pchisq(g1,v) # Significance value of skewness
P1c <- 1 - pchisq(g1c,v) # Significance value of skewness corrected for small sample
g2 <- (b2p-(p*(p+2))) / (sqrt((8*p*(p+2))/n)) # Kurtosis test statistic (approximates to a unit-normal distribution)
P2 <- 1-pnorm(abs(g2)) # Significance value of kurtosis
stats.Hs <- (P1 < alpha);
stats.Ps <- P1;
stats.Ms <- g1;
stats.CVs <- qchisq(1-alpha,v);
stats.Hsc <- (P1c < alpha);
stats.Psc <- P1c;
stats.Msc <- g1c;
stats.Hk <- (P2 < alpha);
stats.Pk <- P2;
stats.Mk <- g2;
stats.CVk <- qnorm(1-alpha,0,1);
cat("Hs = "); print(stats.Hs)
cat("Ps = "); print(stats.Ps)
cat("Ms = "); print(stats.Ms)
cat("CVs = "); print(stats.CVs)
cat("Hsc = "); print(stats.Hsc)
cat("Psc = "); print(stats.Psc)
cat("Msc = "); print(stats.Msc)
cat("Hk = "); print(stats.Hk)
cat("Mk = "); print(stats.Mk)
cat("CVk = "); print(stats.CVk)
}which may be copied into R. When a data matrix, x, is then inputted the program may be run by typing:
mardiatest(x,0.05)