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
|
#
# Wallraff procedure for comparing angular distances
#
# Allows to compare the deviation from an angle of interest
# between several data sets. If the angle of interest is
# the mean direction, then it becomes a comparison of
# angular dispersion around the mean.
#
# In essence, it is a rank-based test (Wilcoxon-Mann-Withney
# or Kruskall-Wallis) on the angular distances from the angle
# of interest.
#
# (c) Copyright 2011 Jean-Olivier Irisson
# GNU General Public License v3
#
#------------------------------------------------------------
# added drop=TRUE 20131106 Claudio
# Generic function
wallraff.test <- function(x, ...) {
UseMethod("wallraff.test", x)
}
# Default method, for an angle vector and a grouping vector
wallraff.test.default <- function(x, group, ref=NULL, ...) {
# 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)")
}
# make sure group is a factor
# if not, force it to keep the order in the original vector
if (!is.factor(group)) {
group <- factor(group, levels=unique(group))
if (!is.null(ref)) {
warning("\"group\" was converted into a factor.\n The levels were kept in the order of the original vector:\n ", paste(levels(group), collapse=", "), "\n Please make sure the elements of \"ref\" match this order")
}
}
# convert data 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")
}
if (is.null(ref)) {
# if no reference angle is provided, use the mean angle
xL <- split(x, group)
ref <- sapply(xL, MeanCircularRad)
} else {
# when a reference angle is provided, if it is not of class circular, cast it as circular assuming that it is in the same reference as the original data
if (!is.circular(ref)) {
ref <- circular(ref, type=dc$type, units=dc$units, template=dc$template, modulo=dc$modulo, zero=dc$zero, rotation=dc$rotation)
}
# now that ref necessarily circular, convert it to the radians/trigonometric case
ref <- conversion.circular(ref, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(ref, "class") <- attr(ref, "circularp") <- NULL
}
# compute concentration parameters and check assumptions
result <- WallraffTestRad(x, group, ref)
result$data.name <- data.name
return(result)
}
# Method for a list
wallraff.test.list <- function(x, ref=NULL, ...) {
# 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)
group <- factor(group, levels=unique(group))
x <- do.call("c", x)
# NB: unlist() removes the circular attributes here
# call default method
result <- wallraff.test.default(x, group, ref)
result$data.name <- data.name
return(result)
}
# Method for a formula
wallraff.test.formula <- function(formula, data, ref=NULL, ...) {
# 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 <- wallraff.test.default(d[,1], d[,2], ref)
result$data.name <- data.name
return(result)
}
# Computation in the usual trigonometric space
WallraffTestRad <- function(x, group, ref) {
# consolidate data
if (length(ref) < nlevels(group)) {
ref = rep(ref, nlevels(group))
}
d <- matrix(c(x, ref=ref[as.numeric(group)]), ncol=2)
# compute angular distances = ranges
dists <- apply(d, 1, function(X) RangeCircularRad(X[1:2], test=FALSE) )
result <- kruskal.test(dists, group)
# NB: kruskal.test with 2 groups is equivalent to wilcox.test with exact=FALSE and correct=FALSE
result$method <- "Wallraff rank sum test of angular distance"
return(result)
}
|
