[R] G-test

[Michael Camann]

Don't know whether anyone ever answered the question regarding G-tests,
but here's a version I cobbled together for another project in case
there's still any interest in this.  This is the G-test ala Sokal and
Rohlf pp 737-738 in my old edition (not in hand, so I don't know the ed.
number).  The first function actually calculates the G statistic and the
second one performs a G-test on a 2x2 table, taking advantage of
print.htest to format the output. They're a bit kludgy, but they seem to
work.  Hope this is of some use....

--Mike C.


# Calculate a G statistic for a 2 x 2 table.
# Continuity corrections are "none", "yates", and "williams", with "none"
# as the default.
#
g.statistic <- function(x, correct="none")	# x is a 2 x 2 matrix
{	
	dname <- deparse(substitute(x))		# the input data
	if (is.data.frame(x)) x <- as.matrix(x)	# coerce data type?

	calc <- function(y) return(y*log(y))	# calculates x ln x

	# reference Sokal and Rohlf pp. 737-738
	# calculate margin sums
	m1 <- apply(x,2,sum)			# column margin sums
	m2 <- apply(x,1,sum)			# row margin sums
	n <- sum(x)				# combined margin sum

	if(correct=="yates")			# Yate's correction?
	{
		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
		}
	}

	# quantity 1
	q1 <- calc(x[1,1]) + calc(x[1,2]) + calc(x[2,1]) + calc(x[2,2])

	# quantity 2
	q2 <- calc(m2[1]) + calc(m2[2]) + calc(m1[1]) + calc(m1[2])
	
	# quantity 3
	q3 <- calc(n)

	# compute G
	G <- 2*(q1 - q2 + q3)

	if (correct=="williams")		# Williams's correction?
	{
		q <- 1+((((n/(m2[1]))+(n/(m2[2]))-1)*((n/(m1[1]))+(n/(m1[2]))-1))/(6*n))
		G <- G/q
	}

	return(G)
}


# this function performs a G test for independence on a 2 x 2 contingency
# table
#
g.test <- function(x, correct="none")		# x is a 2 x 2 matrix
{
	DNAME <- deparse(substitute(x))		# the input data name
	if (is.data.frame(x)) x <- as.matrix(x)	# coerce data type?

	if(dim(x)[1]!=2 || dim(x)[2]!=2)	# insure a 2 x 2 matrix
		stop("x must be a 2 x 2 matrix")

	STATISTIC <- g.statistic(x,correct)	# calc G statistic
	names(STATISTIC) <- "G-statistic"
	PARAMETER <- 1				# X-squared df = 1
	names(PARAMETER) <- "X-squared df"
	PVAL <- 1-pchisq(STATISTIC,df=PARAMETER)# calculate pvalue
	names(PVAL) <- "p.value"

	if(correct=="none") 
		METHOD <- "G-test for independence without continuity 
correction"
	else if(correct=="yates") 
		METHOD <- "G-test for independence with Yate's continuity
correction"
	else if(correct=="williams") 
		METHOD <- "G-test for independence with Williams'
continuity correction"


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


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Michael A. Camann                                  Voice: 707-826-3676
Associate Professor of Zoology                       Fax: 707-826-3201
Institute for Forest Canopy Research     Email: mac24 at axe.humboldt.edu
Department of Biology                            ifcr at axe.humboldt.edu
Humboldt State University           
Arcata, CA 95521

                 URL:http://www.humboldt.edu/~mac24/
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._