R FUN Alexander-Govern

Aus Metheval - Tipps+Tools - Wiki
Wechseln zu: Navigation, Suche

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)