A permutation test appears to work well.
Ordinarily one would think of regressing $B$ and $C$ against $A$. However, since all three are measurements of subjects, all three likely have random components with comparable amounts of dispersion, whereas the regression would assume $A$ has no random component at all. So let's examine this from first principles.
Since you are using correlation to assess relationships among $A$, $B$, and $C$, and since they are purportedly measuring one attribute $X$ of the subjects, it seems fair to suppose those relationships are linear. This suggests modeling the data in terms of the "latent" variable $X$, independent variations $G$, $D$, and $E$, and parameters $(\lambda_0, \ldots, \nu)$, with
$$(A,B,C) = (\lambda_0 + \lambda X + G, \mu_0 + \mu X + D, \nu_0 + \nu X + E).$$
We need to work out the correlations. That starts with the variances and covariances. Writing $\text{Var}(X) = \sigma^2$, $\text{Var}(G) = \sigma_G^2$, etc,
$$\text{Cov}(A,B) = \lambda \mu \sigma^2$$
$$\text{Cov}(A,C) = \lambda \nu \sigma^2$$
$$\text{Var}(A) = \lambda^2 \sigma^2 + \sigma^2_G$$
$$\text{Var}(B) = \mu^2 \sigma^2 + \sigma^2_D$$
$$\text{Var}(C) = \nu^2 \sigma^2 + \sigma^2_E.$$
Therefore the correlations are
$$\eqalign{
\rho(A,B) &= \frac{\text{Cov}(A,B)}{\sqrt{\text{Var}(A)\text{Var}(B)}} \\
&= \frac{\lambda \mu \sigma^2}{\sqrt{\left(\lambda^2 \sigma^2 + \sigma^2_G\right)\left(\mu^2 \sigma^2 + \sigma^2_D\right)}} \\
&= \frac{1}{\sqrt{\left(1 + \sigma^2_G/(\lambda^2\sigma^2)\right)\left(1 + \sigma^2_D/(\mu^2 \sigma^2)\right)}},
}$$
etc. The equality $\rho(A,B)=\rho(A,C)$ therefore is equivalent to
$$\frac{\sigma_D^2}{\mu^2} = \frac{\sigma_E^2}{\nu^2},\ \text{Sgn}(\mu) = \text{Sgn}(\nu).$$
It is unclear how to test this without making distributional assumptions. However, we can take a more general approach by observing that if $D/\mu$ and $E/\nu$ were to have identical distributions, then $B$ and $C$ would have the same distribution up to scale and location, since $(B-\mu_0)/\mu = D/\mu$ would have the same distribution as $(C-\nu_0)/\nu = E/\nu$, and this would imply $\rho(A,B) = \rho(A,C)$. Why not adopt this as the null hypothesis?
If that were the case, then the standardized versions of $B$ and $C$ would be interchangeable conditional on $X$. This would allow us to estimate the null distribution of any statistic that measures the discrepancy between $\rho(A,B)$ and $\rho(A,C)$, such as the absolute value of their difference,
$$T(A,B,C) = |\rho(A,B) - \rho(A,C)|.$$
We may therefore conduct a permutation test simply by standardizing the observations of $B$ and $C$, interchanging their values for a random subset of the subjects, recomputing their correlations, and finding $T$. The observed value of $T$ will be referred to this permutation distribution to determine its significance or lack thereof. In particular, we can compute a permutation p-value as chance that such a random permutation gives a value of $T$ exceeding the observed value.
Although with $n$ subjects there are $2^n$ possible such permutations, which can be huge, we can estimate the permutation distribution by sampling from these permutations. Even a sample of around $500$ or so will work decently.
Here is the R code that results from these considerations. Its argument is my.data, a three-column matrix or data frame in which the first column represents observations of $A$ and the other two columns, observations of $B$ and $C$ (in either order). The second argument n.iter specifies how many permutations to perform. It returns a list consisting of the test statistic $T$, its estimated p-value, and the full array of n.iter simulation results.
cortest <- function(my.data, n.iter=500) {
stat <- function(a, b, c) abs(cor(a, b) - cor(a, c)) # Test statistic function
d <- scale(my.data) # Standardize all variables
n <- dim(d)[1] # Find the number of observations
x <- d[, 1] # Reference the first standardized variable
s <- stat(x, d[,2], d[,3]) # Compute the test statistic
sim <- replicate(n.iter, { # Estimate the permutation distribution:
i <- runif(n) < 1/2 # Choose which values to swap
y <- d[, 2]; y[i] <- d[i, 3] # Swap them
z <- d[, 3]; z[i] <- d[i, 2]
stat(x, y, z) # Compute the test statistic
})
p.value <- mean(c(s, sim) >= s) # Find the estimated p-value
return(list(stat=s, p.value=p.value, sim=sim))
}
To illustrate, let's make some data that exhibit a slight difference in correlations. The inputs are the (true) values of the latent variable $X$ and the parameters $\pi_0 = (\lambda_0, \mu_0, \nu_0)$, $\pi_1 = (\lambda,\mu,\nu)$, and $\sigma = (|\sigma_G/\lambda|, |\sigma_D/\mu|, |\sigma_E/\nu|)$.
make.data <- function(x, pi.0, pi.1, sigma) {
n <- length(x)
data.frame(A=pi.0[1] + pi.1[1]*(x + rnorm(n, sd=sigma[1])),
B=pi.0[2] + pi.1[2]*(x + rnorm(n, sd=sigma[2])),
C=pi.0[3] + pi.1[3]*(x + rnorm(n, sd=sigma[3])))
}
set.seed(17)
x <- rnorm(60)
pi.0 <- rnorm(3)
pi.1 <- abs(rnorm(3))
sigma <- c(1/10, 1/5, 1/3)
my.data <- make.data(x, pi.0, pi.1, sigma)
cor(my.data)
The correlation coefficient estimates shown in the output are $\hat\rho(A,B)=0.9742$ and $\hat\rho(A,C)=0.9449$. Let's run the permutation test and display its results graphically:
sim <- cortest(my.data, 1e4)
par(mfrow=c(1,3))
plot(my.data$A, my.data$B, xlab="A", ylab="B")
plot(my.data$A, my.data$C, , xlab="A", ylab="C")
hist(sim$sim, xlim=range(c(sim$sim, sim$stat)),
xlab="Correlation difference",
main="Simulation Histogram",
sub=paste0("p = ", round(sim$p.value, 3)))
abline(v=sim$stat, lwd=2, col="Red")
The vertical red line is the test statistic $T$ (the absolute difference of the two correlations). The p-value of $0.038$ is the proportion of the histogram to the right of that line. Its small value says this sample of $n=60$ subjects provides evidence of a true underlying difference in the correlations (and indeed there is a slight one, by design, since $|\sigma_D/\mu| = 1/5 \ne 1/3 = |\sigma_E/\nu|$).
We should, at the very least, check that the p-values have a uniform distribution under the null hypothesis. Then we should see how they are distributed under various alternatives. This would be an extensive study, but here are three examples to show how it would be carried out. (Running this will take a minute.)
#
# Null distribution.
#
x <- rnorm(60)
p <- replicate(1e3, {
my.data <- make.data(x, pi.0, pi.1, c(0.1, 0.1, 0.1))
cortest(my.data, 99)$p.value #$
})
hist(p, main="Null distribution")
#
# Some alternatives.
#
p <- replicate(1e3, {
pi.0 <- rnorm(3)
pi.1 <- abs(rnorm(3))
my.data <- make.data(x, pi.0, pi.1, c(1/10, 1/5, 1/3))
cortest(my.data, 99)$p.value #$
})
hist(p, main="Large correlations")
p <- replicate(1e3, {
pi.0 <- rnorm(3)
pi.1 <- abs(rnorm(3))
my.data <- make.data(x, pi.0, pi.1, c(1, 1, 5))
cortest(my.data, 99)$p.value #$
})
hist(p, main="Small correlations")
The test is behaving nicely, at least for a sample size of $60$: the p-values are approximately uniform under (an instance of) the null distribution and they are concentrated at small values for two different deviations from the null.
