# R/cov.llrtest.R In fhernanb/model: This package presents useful functions for modeling regresion.

#### Documented in cov.llrtest

#' Log-likelihood ratio test for the equality of covariance matrices.
#'
#' \code{cov.llrtest} function performs the Log-likelihood ratio
#' test for the equality of at least two covariance matrices, the null
#' hypothesis is \eqn{H0: \Sigma_1=\Sigma_2=\ldots=\Sigma_k}.
#'
#' @param x Is the data set.
#' @param group Is a vector indicating the groups of data set.
#' @param alpha Is the significance level.
#'
#' @examples
#'
#' ################# Test with equal variance
#' mean1 <- c(-2, 4)                   # Mean 1
#' mean2 <- c( 2, 2)                   # Mean 2
#' Sigma <- matrix(c(1, 0.5,
#'                   0.5, 1), ncol=2)  # Covariance matrix
#' library(MASS)
#' n <- 500
#' y1 <- mvrnorm(n, mu=mean1, Sigma=Sigma)
#' y2 <- mvrnorm(n, mu=mean2, Sigma=Sigma)
#'
#' dt <- data.frame(rbind(y1, y2),
#'                  group=rep(c('Pop 1', 'Pop 2'), each=n))
#' with(dt, plot(X1, X2, col=group))
#'
#' cov.llrtest(x=dt[, -3], group=dt$group) #' #' #' ################# Test with unequal variance #' mean1 <- c(-2, 4) # Mean 1 #' mean2 <- c( 2, 2) # Mean 2 #' Sigma1 <- matrix(c(1, 0.5, #' 0.5, 1), ncol=2) # Covariance matrix 1 #' Sigma2 <- matrix(c(1.5, 1, #' 1, 2), ncol=2) # Covariance matrix 2 #' #' n <- 500 #' y1 <- mvrnorm(n, mu=mean1, Sigma=Sigma1) #' y2 <- mvrnorm(n, mu=mean2, Sigma=Sigma2) #' #' dt <- data.frame(rbind(y1, y2), #' group=rep(c('Pop 1', 'Pop 2'), each=n)) #' with(dt, plot(X1, X2, col=group)) #' #' cov.llrtest(x=dt[, -3], group=dt$group)
#'
#' ################# Test with unequal variance and 3 groups
#' mean1 <- c(-2, 4)                    # Mean 1
#' mean2 <- c( 2, 2)                    # Mean 2
#' mean3 <- c( 0, 8)                    # Mean 3
#' Sigma1 <- matrix(c(1, 0.5,
#'                    0.5, 1), ncol=2)  # Covariance matrix 1
#' Sigma2 <- matrix(c(1.5, 1,
#'                    1, 2), ncol=2)    # Covariance matrix 2
#' Sigma3 <- matrix(c(1.5, -1,
#'                    -1, 2), ncol=2)   # Covariance matrix 3
#'
#' n <- 500
#' y1 <- mvrnorm(n, mu=mean1, Sigma=Sigma1)
#' y2 <- mvrnorm(n, mu=mean2, Sigma=Sigma2)
#' y3 <- mvrnorm(n, mu=mean3, Sigma=Sigma3)
#'
#' dt <- data.frame(rbind(y1, y2, y3),
#'                  group=rep(c('Pop 1', 'Pop 2', 'Pop 3'),
#'                            each=n))
#' with(dt, plot(X1, X2, col=group))
#'
#' cov.llrtest(x=dt[, -3], group=dt\$group)
#'
#' @details This test is the multivariate generalization of Bartlett's
#' test of homogeneity of variances. The function was taken and adapted from
#' \url{http://people.stat.sc.edu/hansont/stat730/paketo-libre.pdf},
#' visit the url to obtain more details.
#'
#' @seealso \code{\link{cov.Mtest}} is another test for equality of covariance matrices.
#'
#' @return A list with the following components:
#' \item{test}{The value of the statistic.}
#' \item{degree}{The degrees of the freedom.}
#' \item{critical}{The critical value for the test.}
#' \item{p.value}{The p-value of the test.}
#'
#' @export
#'
cov.llrtest <- function(x, group, alpha = 0.05) {
# x is the data set
# group is a numeric vector indicating the groups of
# the data set a is the significance level, set to 0.05 by default
x <- as.matrix(x)
p <- ncol(x)  # dimension of the data set
n <- nrow(x)  # total sample size
group <- as.factor(group)  # to ensure a group factor
k <- nlevels(group) # number of groups
nu <- rep(0, k)  # the sample size of each group will be stored later
pame <- rep(0, k)
# the next 2 'for' functions separate the k groups and extract the
# covariance matrix of each group the way is not the best but it works
nu <- as.vector(table(group))
mat <- mat1 <- array(dim = c(p, p, k))
# the next 3 lines create the pooled covariance matrix and calculate
# the covariance matrix of each group
j <- 0 # auxiliar variable to count
for (i in levels(group)) {
j <- j + 1
mat[, , j] <- ((nu[j] - 1)/nu[j]) * cov(x[group == i, ])
mat1[, , j] <- (nu[j] - 1) * cov(x[group == i, ])
}
Sp <- apply(mat1, 1:2, sum)/n
for (i in 1:k) pame[i] <- det(solve(mat[, , i]) %*% Sp)
test <- sum(nu * log(pame))  # test statistic
df <- 0.5 * p * (p + 1) * (k - 1)  # degrees of freedom of the asymptotic chi-square
p.value <- 1 - pchisq(test, df)  # p-value of the test statistic
crit <- qchisq(1 - alpha, df)  # critical value of the chi-square distribution
list(test = test, degrees = df, critical = crit, p.value = p.value)
}

fhernanb/model documentation built on Sept. 16, 2017, 11 a.m.