Are two Pearson correlation coefficients different?

Question

I am aware of this question here, which mine was listed as a duplicate of, but it does not fully answer my question. It did however help me progress a little further so thanks, I couldn't find it before. The online calculator on the above answer also disagrees (though marginally) with the vassarstats.net calculator so I think that further backs up my reasoning for not using a black box.

So I'll re-explain my problem and hope that this gets opened to answers:

I have two Pearson correlation coefficients which I would like to compare. Each comes from 2 sets of 40 genetically identical lines and are correlations between male and female trait values.

I have managed to do $z$-transformations on my correlation coefficients in R using the function atanh() and replicated that with a home-made function RtoZ <- function (r) 0.5*log((1+r)/(1-r)).

The problem is what to do next: how do I actually test, in R, whether the two correlations are different?

Answer 1

Just in case that someone (else) has to perform a comparison of correlation coefficients on multiple pairs of variables, here's a function based on rg255's helpful reply to copy:

cor.diff.test = function(x1, x2, y1, y2, method="pearson") {
  cor1 = cor.test(x1, x2, method=method)
  cor2 = cor.test(y1, y2, method=method)

  r1 = cor1$estimate
  r2 = cor2$estimate
  n1 = sum(complete.cases(x1, x2))
  n2 = sum(complete.cases(y1, y2))
  fisher = ((0.5*log((1+r1)/(1-r1)))-(0.5*log((1+r2)/(1-r2))))/((1/(n1-3))+(1/(n2-3)))^0.5

  p.value = (2*(1-pnorm(abs(fisher))))

  result= list(
    "cor1" = list(
      "estimate" = as.numeric(cor1$estimate),
      "p.value" = cor1$p.value,
      "n" = n1
    ),
    "cor2" = list(
      "estimate" = as.numeric(cor2$estimate),
      "p.value" = cor2$p.value,
      "n" = n2
    ),
    "p.value.twosided" = as.numeric(p.value),
    "p.value.onesided" = as.numeric(p.value) / 2
  )
  cat(paste(sep="",
          "cor1: r=", format(result$cor1$estimate, digits=3), ", p=", format(result$cor1$p.value, digits=3), ", n=", result$cor1$n, "\n",
          "cor2: r=", format(result$cor2$estimate, digits=3), ", p=", format(result$cor2$p.value, digits=3), ", n=", result$cor2$n, "\n",
          "diffence: p(one-sided)=", format(result$p.value.onesided, digits=3), ", p(two-sided)=", format(result$p.value.twosided, digits=3), "\n"
  ))
  return(result);
}

Answer 2

Once the Fisher's z transformations are done it is just a case of obtaining p-values

# Correlations    
cor.test (df1$a, df1$b, method = "p")
cor.test (df2$a, df2$b, method = "p")

# function to do fisher transformations 
fisher.z<- function (r1,r2,n1,n2) ((0.5*log((1+r1)/(1-r1)))-(0.5*log((1+r2)/(1-r2))))/((1/(n1-3))+(1/(n2-3)))^0.5

# or this (either version will suffice) 
fisher.z<- function (r1,r2,n1,n2) (atanh(r1) - atanh(r2)) / ((1/(n1-3))+(1/(n2-3)))^0.5

#input n and r from correlations manually (two tailed test)
2*(1-pnorm(abs(fisher.z(r1= ,r2= ,n1= ,n2= ))))

See the final four slides of this presentation and the pnorm() function in r.

Answer 3

In case people are still looking for an easy way to compare two $r$ Pearson correlations. There is a function called paired.r in the package psych for R exactly for that.

Usage:

paired.r(xy, xz, yz=NULL, n, n2=NULL,twotailed=TRUE)

Or as a simple example: For r1 and r2 and a sample size of n just do:

paired.r(r1,r2,n=n)