Generate a random variable with a defined correlation to an existing variable

For a simulation study I have to generate random variables that show a prefined (population) correlation to an existing variable $Y$.

I looked into the R packages copula and CDVine which can produce random multivariate distributions with a given dependency structure. It is, however, not possible to fix one of the resulting variables to an existing variable.

Any ideas and links to existing functions are appreciated!

Conclusion: Two valid answers came up, with different solutions:

1. An R script by caracal, which calculates a random variable with an exact (sample) correlation to a predefined variable
2. An R function I found myself, which calculates a random variable with a defined population correlation to a predefined variable
-
See this related question stats.stackexchange.com/questions/13382/… which directly addresses your question (at least the theory side of it). – Macro Aug 31 '11 at 11:59
The following Q is also strongly related & will be of interest: How to generate correlated random numbers (given means variances and degree of correlation). – gung Oct 16 '12 at 16:07

Here's another one: for vectors with mean 0, their correlation equals the cosine of their angle. So one way to find a vector $x$ with exactly the desired correlation $r$, corresponding to an angle $\theta$:

a) get fixed vector $x_1$ and a random vector $x_2$ b) center both vectors (mean 0), giving vectors $\dot{x}_{1}$, $\dot{x}_{2}$ c) make $\dot{x}_{2}$ orthogonal to $\dot{x}_{1}$ (projection onto orthogonal subspace), giving $\dot{x}_{2}^{\perp}$ d) scale $\dot{x}_{1}$ and $\dot{x}_{2}^{\perp}$ to length 1, giving $\bar{x}_{1}$ and $\bar{x}_{2}^{\perp}$ e) $\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$ is the vector whose angle to $\bar{x}_{1}$ is $\theta$, and whose correlation with $\bar{x}_{1}$ thus is $r$. This is also the correlation to $x_1$ since linear transformations leave the correlation unchanged.

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho


For the orthogonal projection $P$, I used the $QR$-decomposition to improve numerical stability, since then simply $P = Q Q'$.

-
 I was trying to rewrite the code into SPSS syntax. I stumble over your QR decomposition which returns 20x1 column. In SPSS I have Gram-Schmidt orthonormalization (which is also a QR decomposition) but unable to replicate your resultant Q column. Can you chew over your QR action to me please. Or indicate some work-around to get the projection. Thanks. – ttnphns Aug 31 '11 at 21:28 @caracal, P <- X %*% solve(t(X) %*% X) %*% t(X) doesn't produce r=0.6, so that's not the work-around. I'm still confused. (I'd be happy to mimic your expression Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE])) in SPSS but don't know how.) – ttnphns Sep 3 '11 at 9:29 @ttnphns Sorry for the confusion, my comment was for the general case. Applying it to situation in the example: Getting the projection matrix via QR-decomposition is just for numerical stability. You can get the projection matrix as $P=X(X'X)^{-1} X'$ if the subspace is spanned by the columns of matrix $X$. In R, you can here write Xctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1]) because the subspace is spanned by the first column of Xctr. The matrix for the projection onto the orthogonal complement is then I-P. – caracal Sep 4 '11 at 10:59 That's fine now. Thanks. Very nice solution. – ttnphns Sep 4 '11 at 13:24

Let $X$ be your fixed variable and you want to generate $Y$ variable that correlates with $X$ by amount $r$. If $X$ is standardized then $Y= r\cdot X+E$, where $E$ is random variable from normal distribution having mean 0 and $sd=\sqrt{1-r^2}$.

Now, if you want to attain the correlation (almost) exactly $r$, you need to provide that $E$ has zero correlation with $X$. This tightening it to zero can be reached by modifying $E$ iteratively. Unfortunately I can't help with R but in case you use SPSS you can find macro FITVAR on my web-page (See "Fit covariates" collection) that will do all the task from start to end for you for any number of variables at once; for example you could train $Y$ to custom correlations with not one but with several $X$s.

-
Thanks for your answer. That is an empirical/ iterative solution I was thinking about as well. For my simulations, however, I need a more analytical solution without a costly fitting procedure. Fortunately, I just found a solution which I will post shortly ... – Felix S Aug 31 '11 at 11:56
This works for generating bivariate normals but does not work for an arbitrary distribution (or any non-'additive' distribution) – Macro Aug 31 '11 at 11:57
This does work for all kinds of distributions. However, due to the rescaling, $X$ no longer has the original distribution. Also, $Y$ will be some mix of the rescaled distribution of $X$ and an independent normal distribution. Regardless, the true correlation between $X$ and $Y$ will be $r$. – Wolfgang Aug 31 '11 at 12:14
@Felix, it is really worth posting the solution that you've found. Post it as an answer to your question. – ttnphns Aug 31 '11 at 13:37
@Wolfgang, that's certainly true - I assumed ttnphns was talking about generating bivariate normals since $E$ was specified as normal (which isn't necessary, in general). I was only pointing out, as you did, that the resulting marginal distribution for $X$ will not be the same as for $Y$, except (under the specification above) in the special case where $X,Y,\sim N(0,1)$. In general it will be the convolution of the densities of $r X$ and $E$. – Macro Aug 31 '11 at 14:12

Here's another computational approach (the solution is adapted from a forum post by Enrico Schumann). According to Wolfgang (see comments), this is computationally identical to the solution proposed by ttnphns.

In contrast to caracal's solution it does not produce a sample with the exact correlation of $\rho$, but two vectors whose population correlation is equal to $\rho$.

Following function can compute a bivariate sample distribution drawn from a population with a given $\rho$. It either computes two random variables, or it takes one existing variable (passed as parameter x) and creates a second variable with the desired correlation:

# returns a data frame of two variables which correlate with a population correlation of rho
# If desired, one of both variables can be fixed to an existing variable by specifying x
getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) {
if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)}
if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!")

C <- matrix(rho, nrow = 2, ncol = 2)
diag(C) <- 1

C <- chol(C)

X2 <- mar.fun(n)
X <- cbind(X1,X2)

# induce correlation (does not change X1)
df <- X %*% C

## if desired: check results
#all.equal(X1,X[,1])
#cor(X)

return(df)
}


The function can also use non-normal marginal distributions by adjusting parameter mar.fun. Note, however, that fixing one variable only seems to work with a normally distributed variable x! (which might relate to Macro's comment).

Also note that the "small correction factor" from the original post was removed as it seems to bias the resulting correlations, at least in the case of Gaussian distributions and Pearson correlations (also see comments).

-
It seems this is only an approximate solution, i.e., the empirical correlation is not exactly equal to $\rho$. Or am I missing something? – caracal Aug 31 '11 at 21:18
It is easy to show that, except for that "small correction to rho" (whose purpose in this context eludes me), this is exactly the same as what ttnphns suggested earlier. The method is simply based on the Choleski decomposition of the correlation matrix to obtain the desired transformation matrix. See, for example: en.wikipedia.org/wiki/…. And yes, this will only give you two vectors whose population correlation is equal to rho. – Wolfgang Sep 1 '11 at 8:15
The "small correction to rho" was in the original post and is described here. Actually, I don't really understand it; but an investigation of 50000 simulated correlations with rho = .3 shows that without the "small correction" an average of r's of .299 is produced, while with the correction an average of .312 (which is the value of the corrected rho) is produced. Therefore I removed that part from the function. – Felix S Sep 2 '11 at 6:53

Just create a random vector and sort until you get desired r.

-
In what situations would this be preferable to the above solutions? – Andy W Oct 5 '11 at 12:23
A situation where a user wants a simple answer. I read a similar question on the r forum, and its the answer that was given. – Adam Oct 5 '11 at 14:22
Unfortunately this solution is not only computationally inefficient and approximate, it will often fail altogether unless some analysis is first applied to determine an appropriate distribution for the "random vector.". I think there is merit to the underlying idea of just throwing some random numbers at the problem and randomly permuting them (not "sorting" them!) until an approximate $r$ is attained (because this is quick and easy to program), but that idea is not clearly expressed in this short reply. – whuber Oct 5 '11 at 15:12