Skip to content

Instantly share code, notes, and snippets.

@friendly
Created February 10, 2016 15:12
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save friendly/749b5a69a067e02b87dd to your computer and use it in GitHub Desktop.
Box's M-test for Homogeneity of Covariance Matrices
# Box's M-test for Homogeneity of Covariance Matrices
boxM <-
function(mod, ...) UseMethod("boxM")
boxM.default <- function(Y, group)
{
dname <- deparse(substitute(Y))
if (!inherits(Y, c("data.frame", "matrix")))
stop(paste(dname, "must be a numeric data.frame or matrix!"))
if (length(group) != nrow(Y))
stop("incompatible dimensions!")
Y <- as.matrix(Y)
group <- as.factor(as.character(group))
p <- ncol(Y)
nlev <- nlevels(group)
lev <- levels(group)
dfs <- tapply(group, group, length) - 1
if (any(dfs < p))
warning("there are one or more levels with less observations than variables!")
mats <- aux <- list()
for(i in 1:nlev) {
mats[[i]] <- cov(Y[group == lev[i], ])
aux[[i]] <- mats[[i]] * dfs[i]
}
names(mats) <- lev
pooled <- Reduce("+", aux) / sum(dfs)
logdet <- log(unlist(lapply(mats, det)))
minus2logM <- sum(dfs) * log(det(pooled)) - sum(logdet * dfs)
sum1 <- sum(1 / dfs)
Co <- (((2 * p^2) + (3 * p) - 1) / (6 * (p + 1) *
(nlev - 1))) * (sum1 - (1 / sum(dfs)))
X2 <- minus2logM * (1 - Co)
dfchi <- (choose(p, 2) + p) * (nlev - 1)
pval <- pchisq(X2, dfchi, lower.tail = FALSE)
means <- aggregate(Y, list(group), mean)
rn <- as.character(means[,1])
means <- means[,-1]
means <- rbind( means, colMeans(Y) )
rownames(means) <- c(rn, "_all_")
logdet <- c(logdet, pooled=log(det(pooled)))
out <- structure(
list(statistic = c("Chi-Sq (approx.)" = X2),
parameter = c(df = dfchi),
p.value = pval,
cov = mats, pooled = pooled, logDet = logdet, means = means,
data.name = dname,
method = " Box's M-test for Homogeneity of Covariance Matrices"
),
class = c("htest", "boxM")
)
return(out)
}
boxM.formula <- function(formula, data, ...)
{
form <- formula
mf <- model.frame(form, data)
if (any(sapply(2:dim(mf)[2], function(j) is.numeric(mf[[j]]))))
stop("Box's M test is not appropriate with quantitative explanatory variables.")
Y <- mf[,1]
if (dim(Y)[2] < 2) stop("There must be two or more response variables.")
if(dim(mf)[2]==2) group <- mf[,2]
else {
if (length(grep("\\+ | \\| | \\^ | \\:",form))>0) stop("Model must be completely crossed formula only.")
group <- interaction(mf[,2:dim(mf)[2]])
}
boxM.default(Y=Y, group=group, ...)
}
boxM.lm <- function(y, ...) {
boxM.formula(formula(y), data=model.frame(y), ...)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment