[R] G-test : log-likelihood ratio test

[Peter L. Hurd]

I've written a g.test() aka log-likelihood ratio test function for my
opwn use.  It's something I've seen requested (and looked to find
myself) on this list a few times.

It has the same basic syntax as chisq.test().
It does both goodness of fit tests and tests of independence.
Yates' and Williams' corrections are implemented.

I've put some examples from Sokal & Rohlf up at
http://www.psych.ualberta.ca/~phurd/cruft

Feedback welcomed,
-P.

# Kludgy log-likelihood test of independence & goodness of fit
# Does Williams' and Yates' correction
#
# G (TOI) calculation from Zar (2000) Biostatistical Analysis 4th ed.
# G (GOF) calculation from Sokal & Rohlf (1995) Biometry 3rd ed.
# q calculation from Sokal & Rohlf (1995) Biometry 3rd ed.
# TOI Yates' correction taken from Mike Camann's 2x2 G-test fn.
# GOF Yates' correction as described in Zar (2000)
# more stuff taken from ctest's chisq.test()
#
# ToDo:
# 1) Beautify
# 2) Add warnings for violations
# 3) Make appropriate corrections happen by default
#
# V2.1 Pete Hurd Sept 20 2001.

g.test <- function(x, y = NULL, correct="none", p = rep(1/length(x), length(x)))
{
  DNAME <- deparse(substitute(x))
  if (is.data.frame(x)) x <- as.matrix(x)
  if (is.matrix(x)) {
    if (min(dim(x)) == 1) 
      x <- as.vector(x)
  }
  if (!is.matrix(x) && !is.null(y)) {
    if (length(x) != length(y)) 
      stop("x and y must have the same length")
    DNAME <- paste(DNAME, "and", deparse(substitute(y)))
    OK <- complete.cases(x, y)
    x <- as.factor(x[OK])
    y <- as.factor(y[OK])
    if ((nlevels(x) < 2) || (nlevels(y) < 2)) 
      stop("x and y must have at least 2 levels")
    x <- table(x, y)
  }
  if (any(x < 0) || any(is.na(x))) 
    stop("all entries of x must be nonnegative and finite")
  if ((n <- sum(x)) == 0) 
    stop("at least one entry of x must be positive")
  #If x is matrix, do test of independence
  if (is.matrix(x)) {
    #this block was the separate g.stat function
    cell.tot <- row.tot <- col.tot <- grand.tot <- 0
    nrows<-nrow(x)
    ncols<-ncol(x)
    if (correct=="yates"){ # Do Yates' correction
      if(dim(x)[1]!=2 || dim(x)[2]!=2) # check for 2x2 matrix
        stop("Yates' correction requires a 2 x 2 matrix")
      if((x[1,1]*x[2,2])-(x[1,2]*x[2,1]) > 0)
        {
          x[1,1] <- x[1,1] - 0.5
          x[2,2] <- x[2,2] - 0.5
          x[1,2] <- x[1,2] + 0.5
          x[2,1] <- x[2,1] + 0.5
        }
      else
        {
          x[1,1] <- x[1,1] + 0.5
          x[2,2] <- x[2,2] + 0.5
          x[1,2] <- x[1,2] - 0.5
          x[2,1] <- x[2,1] - 0.5
        }
    }
    # calculate G (Zar, 2000)
    for (i in 1:nrows){
      for (j in 1:ncols){
        if (x[i,j] != 0) cell.tot <- cell.tot + x[i,j] * log10(x[i,j])
      }
    }
    for (i in 1:nrows){ row.tot <- row.tot + (sum(x[i,])) * log10(sum(x[i,])) }
    for (j in 1:ncols){ col.tot <- col.tot + (sum(x[,j])) * log10(sum(x[,j])) }
    grand.tot <- sum(x)*log10(sum(x))
    total <- cell.tot - row.tot - col.tot + grand.tot
    G <- 4.60517 * total
    q <- 1
    if (correct=="williams"){ # Do Williams' correction
      row.tot <- col.tot <- 0    
      for (i in 1:nrows){ row.tot <- row.tot + 1/(sum(x[i,])) }
      for (j in 1:ncols){ col.tot <- col.tot + 1/(sum(x[,j])) }
      q <- 1+ ((n*row.tot-1)*(n*col.tot-1))/(6*n*(ncols-1)*(nrows-1))
    }
    G <- (G/q)
    # end of old g.stat function
    
    STATISTIC <- G
    PARAMETER <- (nrow(x)-1)*(ncol(x)-1)
    PVAL <- 1-pchisq(STATISTIC,df=PARAMETER)
    if(correct=="none")
      METHOD <- "Log likelihood ratio (G-test) test of independence without correction"
    if(correct=="williams")
      METHOD <- "Log likelihood ratio (G-test) test of independence with Williams' correction"
    if(correct=="yates")
      METHOD <- "Log likelihood ratio (G-test) test of independence with Yates' correction"
  }
  else {
    # x is not a matrix, so we do Goodness of Fit
    METHOD <- "Log likelihood ratio (G-test) goodness of fit test"
    if (length(x) == 1) 
      stop("x must at least have 2 elements")
    if (length(x) != length(p)) 
      stop("x and p must have the same number of elements")
    E <- n * p
    
    if (correct=="yates"){ # Do Yates' correction
      if(length(x)!=2)
        stop("Yates' correction requires 2 data values")
      if ( (x[1]-E[1]) > 0.25) {
        x[1] <- x[1]-0.5
        x[2] <- x[2]+0.5
      }
      else if ( (E[1]-x[1]) > 0.25){
        x[1] <- x[1]+0.5
        x[2] <- x[2]-0.5
      }
    }
    names(E) <- names(x)
    tot <- 0
    for (i in 1:length(x)){
      if (x[i] != 0) tot <- tot + x[i] * log(x[i]/E[i])
    }
    G <- (2*tot)
    if (correct=="williams"){ # Do Williams' correction
      q <- 1+(length(x)+1)/(6*n)
      G <- (G/q)
    }
    STATISTIC <- (G)
    PARAMETER <- length(x) - 1
    PVAL <- pchisq(STATISTIC, PARAMETER, lower = FALSE)
  }
  names(STATISTIC) <- "Log likelihood ratio statistic (G)"
  names(PARAMETER) <- "X-squared df"
  names(PVAL) <- "p.vlaue"
  structure(list(statistic=STATISTIC,parameter=PARAMETER,p.value=PVAL,
                 method=METHOD,data.name=DNAME),class="htest")
}

-- 
Peter L. Hurd                                  Department of Psychology
Assistant Professor                               University of Alberta
phurd at ualberta.ca                                     Edmonton, Alberta
http://www.psych.ualberta.ca/~phurd                      T6G 2E9 Canada


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._