# Bubba's log-likelihood & chi-squared tests of independence
#   for three dimensional data matrices
#
# ToDo:
# 1) Beautify - why isn't output appropriate for class?
# 2) Add documentation, as in: if called with 'method="c"'
#    then chi-squared is used, else defaults to g-test
# 3) Add corrections
#
# V1.0 Pete Hurd 31 July 2004. phurd@ualberta.ca
#
threeD.toi.test <- function(x,method="g"){
  size <- dim(x)
  size.x <- size[1]
  size.y <- size[2]
  size.z <- size[3]

  sr <- apply(x,1,sum)
  sc <- apply(x,2,sum)
  
  n <- sum(x)
  E <- outer(sr,sc, "*")/n

  nz  <- numeric(length=0)
  for (z in 1:size.z){
    nz <- c(nz,sum(x[,,z]))
  }
  E3 <- numeric(length=0)
  for (i in nz){  E3 <- c(E3,as.vector(E)*i/n)}
  dim(E3) <- size

  if (method=="c"){
    val <- sum((x - E3)^2/E3)
  }
  else{ 
    val <- 0
    for (i  in 1:size.x) {
      for (j  in 1:size.y) {
        for ( k in 1:size.z) {
          if (x[i,j,k] != 0){
            val <- val + x[i,j,k] * log(x[i,j,k]/E3[i,j,k])
          }
        }
      }
    }
  }

  STATISTIC <- val
  print(STATISTIC)
  if (method=="c") names(STATISTIC) <- "X-squared"
  else names(STATISTIC) <- "Log-likelihood ratio"

  METHOD <- "Test of Independence on a three-dimensional matrix"
  PARAMETER <- (size.x - 1)*(size.y - 1)*(size.z - 1) 
  PVAL <- pchisq(STATISTIC, PARAMETER, lower = FALSE)
  DNAME <- deparse(substitute(x))

  structure(list(statistic = STATISTIC, parameter = PARAMETER, 
                 p.value = PVAL, method=METHOD, data.name = DNAME, observed = x, 
                 expected = E3, class = "htest"))
  
}