another function to calculate gtest in R

gTest2.R

#http://www.stat.wisc.edu/~st571-1/gtest.R
#see page 2 for data & page 17 for g-test http://www.stat.wisc.edu/~st571-1/06-tables-2.pdf

# Bret Larget
# October 28, 2010
# September 30, 2011


# Use source() to read this file into R to use it.
# Example:
#   # First, change the working directory in R to the location where the file gtest.R resides.
#   # Then, use source() to read in this file.
#   source("gtest.R")
#   # Create a matrix with data for this test.
#   fish = matrix(c(1,49,10,35,37,9),nrow=2,ncol=3)
#   rownames(fish) = c("Eaten","Not Eaten")
#   colnames(fish) = c("Uninfected","Lightly Infected","Highly Infected")
#   # Call the function.
#   g.test(fish)


# Function to carryout G-test for contingency table
# x is a matrix with (non-negative) integer counts

g.test = function(x,print=T) {
  # find the row and column sums of x with apply
  # apply(x,1,f) applies function f to each row (first dimension) of x
  # apply(x,2,f) applies function f to each column (second dimension) of x
  row.sum = apply(x,1,sum)
  col.sum = apply(x,2,sum)
  # n is the sum of all elements of the table
  n = sum(x)
  # x.expected[i,j] = row.sum[i] * col.sum[j] / n
  x.expected = row.sum %o% col.sum / n
  # g = 2 * sum_i sum_j ( x[i,j] * log(x[i,j]/x.expected[i,j] )
  g = 2*sum( x*log(x/x.expected) )
  # degrees of freedom is (#rows-1)*(#columns-1)
  degf = (nrow(x)-1)*(ncol(x)-1)
  # If conditions are right, p.value is from chi^2 distribution (area to right of test statistic)
  p.g = 1 - pchisq(g,degf)
  # return a list with expected counts, test statistic, and p-value
  if ( print ) {
    cat("G-Test for Contingency Tables\n\n")
    cat("Data:\n")
    print(x)
    cat("\n")
    cat("The test statistic is ",signif(g,6),".\n")
    cat("There are ",degf," degrees of freedom.\n")
    cat("The p-value is ",p.g,".\n\n\n")
  }
  return(invisible(list(expected=x.expected,test.stat=g,p.value=p.g)))
}