Solved – Are two Pearson correlation coefficients different

correlationhypothesis testingpearson-rr

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?

Best Answer

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);
}
Related Question