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
|
###############################################################
# #
# Original Splus: Ulric Lund #
# E-mail: ulund@calpoly.edu #
# #
###############################################################
###############################################################
# #
# R port: Claudio Agostinelli <claudio@unive.it> #
# #
# Date: May, 26, 2006 #
# Version: 0.3-1 #
# #
###############################################################
cor.circular <- function(x, y=NULL, test = FALSE) {
if (!is.null(y) & NROW(x)!=NROW(y))
stop("x and y must have the same number of observations")
if (is.null(y) & NCOL(x)<2)
stop("supply both x and y or a matrix-like x")
ncx <- NCOL(x)
ncy <- NCOL(y)
# Handling missing values
if (is.null(y)) {
ok <- complete.cases(x)
x <- x[ok,]
} else {
ok <- complete.cases(x, y)
if (ncx==1) {
x <- x[ok]
} else {
x <- x[ok,]
}
if (ncy==1) {
y <- y[ok]
} else {
y <- y[ok,]
}
}
n <- NROW(x)
if (n==0) {
warning("No observations (at least after removing missing values)")
return(NULL)
}
x <- conversion.circular(x, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(x, "class") <- attr(x, "circularp") <- NULL
if (!is.null(y)) {
y <- conversion.circular(y, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(y, "class") <- attr(y, "circularp") <- NULL
}
if (is.null(y)) {
result <- matrix(1, ncol=ncx, nrow=ncx)
if (test) {
test.stat <- matrix(0, ncol=ncx, nrow=ncx)
p.value <- matrix(0, ncol=ncx, nrow=ncx)
}
for (i in 1:ncx) {
for (j in i:ncx) {
res <- CorCircularRad(x=x[,i], y=x[,j], test=test)
result[i,j] <- result[j,i] <- res[1]
if (test) {
if (i==j) {
test.stat[i,i] <- NA
p.value[i,i] <- NA
} else {
test.stat[i,j] <- test.stat[j,i] <- res[2]
p.value[i,j] <- p.value[j,i] <- res[3]
}
}
}
}
} else {
attributes(x) <- c(attributes(x), list(dim=c(n, ncx)))
attributes(y) <- c(attributes(y), list(dim=c(n, ncy)))
result <- matrix(1, ncol=ncy, nrow=ncx)
if (test) {
test.stat <- matrix(0, ncol=ncy, nrow=ncx)
p.value <- matrix(0, ncol=ncy, nrow=ncx)
}
for (i in 1:ncx) {
for (j in 1:ncy) {
res <- CorCircularRad(x=x[,i], y=y[,j], test=test)
result[i,j] <- res[1]
if (test) {
test.stat[i,j] <- res[2]
p.value[i,j] <- res[3]
}
}
}
}
if (ncx==1 | (!is.null(y) & ncy==1)) {
result <- c(result)
if (test) {
test.stat <- c(test.stat)
p.value <- c(p.value)
}
}
if (test) {
result <- list(cor=result, statistic=test.stat, p.value=p.value)
}
return(result)
}
CorCircularRad <- function(x, y, test=FALSE) {
n <- length(x)
x.bar <- MeanCircularRad(x)
y.bar <- MeanCircularRad(y)
num <- sum(sin(x - x.bar) * sin(y - y.bar))
den <- sqrt(sum(sin(x - x.bar)^2) * sum(sin(y - y.bar)^2))
result <- num/den
if (test) {
l20 <- mean.default(sin(x - x.bar)^2)
l02 <- mean.default(sin(y - y.bar)^2)
l22 <- mean.default((sin(x - x.bar)^2) * (sin(y - y.bar)^2))
test.stat <- sqrt((n * l20 * l02)/l22) * result
p.value <- 2 * (1 - pnorm(abs(test.stat)))
result <- c(result, test.stat, p.value)
}
return(result)
}
|
