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

Question

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

Answer 1

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'$.

Answer 2

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.

Answer 3

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).

Answer 4

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