# R FUN Alexander-Govern

The Alexander-Govern approximation in R (r-project.org)

```## Alexander-Govern approximation test
## (c) 2007 Sven Hartenstein

## This code comes with NO WARRANTY, not even the warranty that it
## correctly computes the test statistic as described by Alexander and
## Govern. I did not find any software computing the Alexander-Govern
## test statistic to compare the results of this function with. I
## assume however that it is correct and some simulations are in line
## with this assumption.

## Note that the function does NOT check for incorrect call, invalid
## parameters or NAs. It rather assumes that there are no missings and
## that the user calls it appropriately.

## For a description and development of the test statistic see
## Alexander, R.A., & Govern, D.M. (1994). A New and Simpler
## Approximation for ANOVA Under Variance Heterogeneity. Journal of
## Educational Statistics, 19, 91-101.

ag.test <- function(y, group) {
## y should be a vector whichs elements are the values of the
## numeric dependent variable. group should be a vector of the same
## length whichs elements (numeric, character or mixed) indicate the
## group each observation belongs to.

## The function returns a list with the following elements:
## - \$approx: the test statistic (chisq distributed)
## - \$df: the degrees of freedom
## - \$p.value: the p value

## some preliminary computations
n <- length(y)
x.levels <- levels(factor(group))
y.se <- w <- y.means <- y.n <- NULL
for (i in x.levels) {
## compute standard errors (formula (1))
y.se[i] <- sqrt(var(y[group==i])/length(y[group==i]))
## we need the group means below
y.means[i] <- mean(y[group==i])
## we need the group sizes as well
y.n[i] <- length(y[group==i])
}
## compute weights (formula (2))
w <- (1/y.se^2)/sum(1/y.se^2)

## compute variance-estimate of the common mean of Y (formula (3))
y.plus <- sum(w * y.means)

## compute one-sample t statistic (formula (4))
t <- (y.means - y.plus)/y.se

## compute the normalizing transformations of t (formula (8))
a <- y.n - 1.5
b <- 48*a^2
c <- (a * log(1 + t^2/(y.n-1)))^(.5)
z <- c + (c^3+3*c)/b - (4*c^7+33*c^5+240*c^3+855*c)/(10*b^2+8*b*c^4+1000*b)

## compute approximation (formula (9))
approx <- sum(z^2)
df <- length(x.levels)-1
p.value <- pchisq(approx, df=df, lower.tail=FALSE)

## return result
return(list(approx=approx, df=df, p.value=p.value))
}

### Example:
## set.seed(1)
## y <- c(rnorm(40, sd=10),
##        rnorm(30, sd=15),
##        rnorm(20, sd=20))
## x <- c(rep(1, times=40),
##        rep("two", times=30),
##        rep(8, times=20))
## ag.test(y, x)
```