1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
#
# 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
#
#------------------------------------------------------------
# added drop=TRUE 20131106 Claudio
# 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
kt <- EqualKappaTestRad(x, group)
# 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 <- WatsonWilliamsTestRad(x, group, kt)
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
means <- tapply(x, group, MeanCircularRad)
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)
}
|
