# Why does the multivariate data generated by a copula in R not exhibit the prespecified correlation?

I am using the package copula in R to generate a bivariate sample. The marginal distributions are binomial with p=0.5 and exponential with rate=1. With prespecified rho=0.5, the data I generated only has correlation 0.35; actually, as I changed the value of rho, I found that the sample correlation is always 70% of the prespecified rho.

The code is shown below:

library(copula)
library(MASS)
library(psych)

rho = 0.5
cop_control <- normalCopula(param = rho, dim = 2)
bi_control <- mvdc(copula=cop_control, margins=c("binom", "exp"),
paramMargins=list(list(size=1, prob=0.5),list(rate=1)))
u_control <- rMvdc(5000, bi_control)
colnames(u_control) <- c("x1", "x2")
pairs.panels(u_control)


Shouldn't the sample from the copula bivariate distribution have the prespecified correlation?

How should I change my code to accurately sample? Thanks!

Three things to note:

1. Pearson correlation is not preserved through strictly monotonic transformation of the margins, including marginal transformations used with copulas, except in special cases. When it's not linear generally the correlation reduces*. However, monotonic correlations including the Spearman and Kendall correlations are preserved (this is a reason why they're heavily used with copulas).

2. You have in there a not-strictly-monotonic transformation (the transformation to binomial). So even the monotonic correlations are no longer preserved.

3. Sample correlation differs from population correlation and your sample size is only 5000. So you might easily see a correlation not too far from 0 differ by say 0.025 or so from its population value in such a sample. Here's the third sample I generated just now (my examples use $$\rho=0.6$$):

x=rnorm(5000); y=.6*x+.8*rnorm(5000) ; cor(x,y)
[1] 0.5743064


It's worth keeping that last one in mind; sometimes people are surprised that at some large sample size there's still sometimes a substantive difference between sample and population values.

In short, in general you cannot specify a population Pearson correlation, transform the variables nonlinearly and still have the same correlation between the new variables that you started with. Numerous posts on site address this basic issue.

If you want a certain population Pearson correlation post-transformation (in spite of not having E(Y|X=x) being linear in x) you can approximate it with a little trial and error. For example if you're getting the final correlation being 70% of what you wanted, across different initial rho values, as a first attempt, start with rho increased by a multiple of 1/0.7 so that 70% of it is roughly where you want to end up; e.g. if you want to end up at 0.5, start at about 0.714. You might get there with fewer attempts using an approximation to the transformed moments perhaps, though the overall time used to get there is probably longer than just simple trial and error unless you're doing it many times. Once you get within the typical variation in sample values from the population values, you're about as close as you're going to get to the desired population correlation without using a larger sample size.

We can see the effects described in 1 and 2 above:

In the first of the three plots above we see 5000 points from a bivariate normal with rho=0.6 (sample r is 0.601 here, and the sample Kendall tau is also given - it's 0.41). I blur the distinction between sample and population values noted in the plot but the variation is small enough that the point is clear either way.

In the second plot (top right), we see an example of the fact that if we apply nonlinear monotonic transformations, the resulting relationship is nonlinear (unless they're carefully contrived otherwise), and the population and sample Pearson correlations change as a result, but $$\tau$$ values do not (similarly with Spearman correlation, quadrant correlation, etc).

In the third plot (bottom left) I apply the marginal transformations it looks like you're using, but if I have that detail wrong somewhere, the underlying point still stands and the illustration still serves as explanation for what you observe and my discussion of it. Both correlations are now changed because the binomial transformation, while monotonic, is no longer strictly monotonic (we transform from a continuous variable to a discrete one, so multiple values end up "in the same place"; these transformations are not invertible).

* but not necessarily always; e.g. it might just be possible to stretch the ends enough in just the right way that correlation sometimes increases despite the nonlinearity. I have no examples of it, nor proof that it's possible, so perhaps the initial $$\rho$$ is actually an upper bound on the correlation of the transformed values. This difficulty in coming up with an example might relate to the absence of tail dependence in the Gaussian copula; I think it would be possible with a copula with sufficient tail dependence.

No, the correlation that you specify is on the latent Gaussian scale. The correlation will typically (?always) be lower for the observed variables with the specified distributions