# R/cor.circular.R In circular: Circular Statistics

#### Documented in cor.circular

```###############################################################
#                                                             #
#       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) {
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) {
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)
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)
}
```

## 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.