# R/watson.williams.test.R In circular: Circular Statistics

#### Documented in watson.williams.testwatson.williams.test.defaultwatson.williams.test.formulawatson.williams.test.list

#
#	Watson-Williams test for homogeneity of means
#
#	Allows to compare mean angles in two or more samples.
#	Equivalent, for angles, of an ANOVA/Kruskal-Wallis test.
#
#	Based on
#		Circular statistics in biology, Batschelet, E (1981)
#		6.2, p. 99
#		Biostatistical analysis, Zar, J H (1999)
#		27.4, p. 634
#		Directional statistics, Mardia, K.V. and Jupp, P.E. (2000)
#		p. 135
#
# (c) Copyright 2010-2011 Jean-Olivier Irisson
#     GNU General Public License, v.3
#
#------------------------------------------------------------

# Generic function
watson.williams.test <- function(x, ...) {
UseMethod("watson.williams.test", x)
}

# Default method, for an angle vector and a grouping vector
watson.williams.test.default <- function(x, group, ...) {

# get data name
data.name <- paste(deparse(substitute(x)), "by", deparse(substitute(group)))

# check arguments
ok <- complete.cases(x, group)
x <- x[ok]
group <- group[ok,drop=TRUE]
if (length(x)==0 | length(table(group)) < 2) {
stop("No observations or no groups (at least after removing missing values)")
}

# convert everything to the radians/trigonometric case
if (is.circular(x)) {
dc <- circularp(x)
x <- conversion.circular(x, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(x, "class") <- attr(x, "circularp") <- NULL
} else {
dc <- list(type="angles", units="radians", template="none", modulo="asis", zero=0, rotation="counter")
}

# compute concentration parameters and check assumptions

# equality of concentration parameters
if (kt\$p.value < 0.05) {
warning("Concentration parameters (", paste(format(kt\$kappa, digits=3), collapse=", ") ,") not equal between groups. The test might not be applicable")
}

# sufficiently large concentration
# Batschelet provides kappa.all > 2 (or equivalently rho.all > 0.75) in the two sample case (but no indication in the multisample one)
# Mardia & Jupp cite Stephens 1972 to justify that kappa.all >= 1 (or equivalently rho.all >= 0.45) in the multisample case
# Zar's adds conditions on minimum sample size to use the smaller thresholds of the concentration parameter -- but there is no discussion of how things change when the sample size smaller but the concentration larger) :
#   N / 2 >= 25 in the two sample case
#   N / k >= 6 in the multisample case
# determine whether we are doing a two or multisample test
if (length(table(group)) == 2) {
kappa.thresh <- 2
} else {
kappa.thresh <- 1
}
if ( kt\$kappa.all < kappa.thresh ) {
warning("Global concentration parameter: ", format(kt\$kappa.all, digits=3)," < ", kappa.thresh, ". The test is probably not applicable")
}

# TODO : also check that distributions conform to Von Mises?

result\$data.name <- data.name

# convert means back in their original units
result\$estimate <-  conversion.circular(circular(result\$estimate), units=dc\$units, type=dc\$type, template=dc\$template, modulo=dc\$modulo, zero=dc\$zero, dc\$rotation)

return(result)
}

# Method for a list
watson.williams.test.list <- function(x, ...) {
# fecth or fill list names
k <- length(x)
if (is.null(names(x))) {
names(x) <- 1:k
}
# get data name
data.name <- paste(names(x), collapse=" and ")

# convert into x and group
ns <- lapply(x, length)
group <- rep(names(x), times=ns)
x <- do.call("c", x)
# NB: unlist() removes the circular attributes here

# call default method
result <- watson.williams.test.default(x, group)
result\$data.name <- data.name

return(result)
}

# Method for a formula
watson.williams.test.formula <- function(formula, data, ...) {
# convert into x and group
d <- model.frame(as.formula(formula), data)

# get data name
data.name <- paste(names(d), collapse=" by ")

# call default method
result <- watson.williams.test.default(d[,1], d[,2])
result\$data.name <- data.name

return(result)
}

# Computation in the usual trigonometric space
WatsonWilliamsTestRad <- function(x, group, kt) {

# number of groups
group <- as.factor(group)
k <- nlevels(group)
# total sample size
n <- length(x)
# sample size per group
ns <- as.numeric(table(group))
# correction factor
g <- 1 + 3 / (8 * kt\$kappa.all)
# sum of resultant vectors lengths
sRi <- sum(kt\$rho * ns)
# total resultant vector length
R <- kt\$rho.all * n

statistic <- g * ((n - k) * (sRi - R)) / ((k - 1) * (n - sRi))

p.value <- pf(statistic, k-1, n-k, lower.tail=FALSE)

# compute estimates of means
names(means) <- paste("mean of", names(means))

# return result
result <- list(
method = "Watson-Williams test for homogeneity of means",
parameter = c(df1=k-1, df2=n-k),
statistic = c(F=statistic),
p.value = p.value,
estimate = means
)
class(result) <- "htest"

return(result)
}

## Try the circular package in your browser

Any scripts or data that you put into this service are public.

circular documentation built on May 1, 2019, 7:57 p.m.