How to calculate p-value from the expected correlation coefficient using randomisation in R?

Question

Here, I want to estimate the p-value using randomization. Tests of correlations involving patterns like Y (a function of X) vs. X involve a shared term and the X term. Here I want to randomize (999 times) the non­shared variable while keeping the shared variable constant (e.g. I randomly shuffled X, not Y). Then I estimate the expected correlation, E(r) as the average of these randomizations. The P­-value associated with this more conservative test required tallying the number of randomizations, rrand,j producing a difference from the expected correlation, |rrand,j - E(r)| that is larger in magnitude than the difference between the observed correlation and the expected correlation, |robs,j - E(r)|

##DATA FRAME
set.seed(111)
library(truncnorm)
x <- rtruncnorm(n = 288,a = 0,b = 10,mean = 5,sd = 2)
v <- rtruncnorm(n = 288,a = 0,b = 10,mean = 5,sd = 2)
y <- ((v/x^2) - (1/x))
sp <- rep(c("A","B","C","D"), each = 72)

df <- data.frame(v,x,y,sp)


##EXPECTED CORRELATION FUNCTION
library(data.table)
setDT(df)
# function to estimate model coefficients
f <- function(x,v) {x.sample <- sample(x, length(x), replace=T)
y.sample <- (v/x.sample^2) - (1/x.sample)
per <- cor(y.sample, x.sample)}


# 999 models for each species
result = rbindlist(
  lapply(1:999, \(i) df[,.(est = f(x,v)), sp][, i:=i])
)

#MEAN EXPECTED CORRELATION
result.mean.exp.cor <-
  result %>% 
  group_by(sp)%>% 
  dplyr::summarise(mean.expected.cor = mean(est))

#TESTING IF R.RAND - E(R) is greater than 0
rrand.exp <- dplyr::left_join(result,result.mean.exp.cor,by = 'sp')
rrand.exp$count <- ifelse(rrand.exp$est > rrand.exp$mean.expected.cor, 1,0)

#TALLYING THE CASES WHERE R.RAND - E(R) is greater than 0
rrand.exp.pval <-
  rrand.exp %>% 
  group_by(sp)%>% 
  dplyr::summarise(pvals = mean(count))

I understand how to find the |rrand,j - E(r)|.I thought that this in itself is the p-value. I don't understand how or why I need to calculate "that is larger in magnitude than the difference between the observed correlation and the expected correlation, |robs,j - E(r)|" Wouldn't this be a single number (for each species)? Am I missing something?

I also don't understand how to calculate |robs,j - E(r)| and compare it to |rrand,j - E(r)|. In this problem, I don't even know where to start because I am comparing 999 values of |robs,j - E(r)| against one of |rrand,j - E(r)|

For all those asking: I'm trying to avoid spurious correlation issues. More reading at https://link.springer.com/article/10.1007/BF00317404 Jackson and Somers(1991).

Thanks everyone. I was making a silly error on how I was thinking about this problem. I figured it out!

Answer 1

Without getting into whether this test problem is generally valid (see comments of the question), I can say this:

The p-value is a probability. But $|r_{rand,j}-E(r)|$ is not a probability, and neither is their mean (note that I use $j$ here for the specific permutation; therefore in my notation $r_{obs}$ has no $j$ index). Generally the p-value is the probability that under null hypothesis ($E(r)$ correlation) the correlation value is as far away (in absolute value) or farther away from the expected value $E(r)$ than what was observed.

In a permutation test, the null hypothesis is operationalised by the permutation, meaning that the probability of interest, i.e., the p-value, is the relative frequency of how often $|r_{rand,j}-E(r)|$ is larger than or equal to $|r_{obs}-E(r)|$.