powdiv<-function(obs, exp, npar, lambda)
#
# Function to calculate Cressie-Read power divergence goodness-of-fit statistic
#
# Input parameters
#
#   obs     - vector of observed counts
#   exp     - vector of expected counts
#   npar    - number of parameters estimated in null model
#   lambda  - power for statistic
#                lambda = 1    : Pearson statistic
#                lambda = 2/3  : Cressie-Read statistic
#                lambda = 0    : likelihood ratio statistic
#                lambda = -1/2 : Freeman-Tukey statistic
#                lambda = -1   : modified likelihood ratio statistic
#                lambda = -2   : Neyman-modified statistic
#
# Output parameters
#
#   stat    - statistic
#   df      - degrees of freedom for statistic
#   pval    - p-value for statistic based on a chi-squared distribution with
#             degrees of freedom equal to df
#
# Note: For all statistics except the Pearson statistic, cells with observed counts
#       equal to zero do not contribute to the test statistic (effectively this
#       means that 0 * log(0) and 0 * Infinity are taken to be zero). This can lead to
#       strange values of the statistic (such as negative values) for lambda less
#       than zero.
#
{
        obs<-as.vector(obs)
        exp<-as.vector(exp)
	if(lambda == 0) {
		stat <- 2 * sum(na.omit(obs * log(obs/exp)))
	}
	else {
		if(lambda == -1) {
			stat <- 2 * sum((exp * log(exp/obs))[
				is.finite(exp * log(exp/obs))])
		}
		else {
			stat <- 2/(lambda * (lambda + 1)) * sum(
				na.omit(obs * ((obs/exp)^lambda - 1
				)))
		}
	}
	df <- length(obs) - npar - 1
	pval <- 1 - pchisq(stat, df)
	return(list(stat = stat, df = df, pval = pval))
}