R/cov.llrtest.R

Defines functions cov.llrtest

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.