I'm attempting to check the correlation coefficient between an angular/circular variable and a linear variable using the circular package, as described in this source: https://digitalcommons.calpoly.edu/cgi/viewcontent.cgi?article=1056&context=statsp. Specifically I'm using the parametric method.
First I created two vectors and coerced one to circular:
theta = c(45,36,57,89,16,39,50,48,41,46)
y = as.circular(theta,units = "degrees",
type = "angles")
x = c(10,8,8,4,6,8,9,9,9,10)
I then copy/pasted the following function from the appendix of the source I linked above:
#####Function#####
cor.circular.lc = function(x,y=NULL,test =
FALSE){
### x vector or matrix of linear data
### y vector or matrix of circular data
### test if test == TRUE then a
### significance test for the correlation
### is computed
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)
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)
}
#### Converting y to radians ####
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)){
z = conversion.circular(x[,2], units =
"radians", zero = 0,
rotation =
"counter",
modulo = "2pi");
attr(z, "class") <- attr(z, "circularp") <- NULL;
r_xs = cor(x[,1],sin(z));
r_xc = cor(x[,1],cos(z));
r_cs = cor(cos(z),sin(z));}else{
#### calculating individual components ####
r_xs = cor(x,sin(y));
r_xc = cor(x,cos(y));
r_cs = cor(cos(y),sin(y));}
#### calculating correlation coeff linear-circular ####
cor.lc = (r_xc^2 + r_xs^2 - 2*(r_xc*r_xs*r_cs))/(1-r_cs^2);
if(test){
f.stat = (.5*(n-3)*cor.lc)/(1-cor.lc);
p.val = pf(f.stat,df1 = 2, df2=
n-3,lower.tail = FALSE);
result = list(cor = cor.lc, statistic =
f.stat, p.value = p.val);
}else{
result = list(cor = cor.lc);
}
return(result);
}
But when I try to run the function with the variables, instead of returning output with three values, it only returns one value, leaving out the F statistic and p-value:
> cor.circular.lc(x,y)
$cor
[1] 0.8858006
From what I can tell, this is being caused by the }else{result = list(cor = cor.lc); line causing only cor to be returned and not f.stat or p.val. The note at the beginning of the function (#test if test == TRUE then a significance test for the correlation is computed) suggests to me that currently test == FALSE for some reason, and when I look at the first few lines of the function, I see function(x,y=NULL,test = FALSE) which further suggests to me that the y is NULL for some reason. Is my logic here sound? If so, how do I fix it? If not, what's going wrong? I need the F statistic and p-value for this test.
(this is my first time asking a question on here, so let me know if there's something I've forgotten to include)
function(x,y=NULL,test = FALSE) means that FALSE is the default value for test. If you want to compute the significance test you should include test = TRUE in the function call:
>cor.circular.lc(x, y, test = TRUE)
$cor
[1] 0.8858006
$statistic
[1] 27.14815
$p.value
[1] 0.0005032954