Create 3 correlated variables with pre-specified correlations I start with a vector of random numbers sampled from a normal distribution:
R<-rnorm(100, mean=0, sd=30)

I would now like to create 3 variables that are correlated with each other with a pre-specified correlation. In addition I would like to have these three variables correlated with R with a pre-specified correlation.
For example A, B, C would have correlation 0.7 with each other, and A, B, C would have correlation 0.6 with R.
I.e. I am looking for the following covariance matrix:
    R   A   B   C
R   1  .6  .6  .6
A  .6   1  .7  .7
B  .6  .7   1  .7
C  .6  .7  .7   1

How can this be done in R?
 A: You could do something like this:
library(mvtnorm)

x = 3.3
sig = matrix(c(30,x,x,x,x,1,.7,.7,x,.7,1,.7,x,.7,.7,1),nrow=4)

X = rmvnorm(100,mean=rep(0,4),sigma=sig,method="svd")
round(cor(X),2)

Y = rmvnorm(10000,mean=rep(0,4),sigma=sig,method="svd")
round(cor(Y),2)

Not the structure of the covariance matrix:
> sig
     [,1] [,2] [,3] [,4]
[1,] 30.0  3.3  3.3  3.3
[2,]  3.3  1.0  0.7  0.7
[3,]  3.3  0.7  1.0  0.7
[4,]  3.3  0.7  0.7  1.0

Now, as you can see with only 100 samples the calculated correlation is:
> X = rmvnorm(100,mean=rep(0,4),sigma=sig,method="svd")
> round(cor(X),2)
     [,1] [,2] [,3] [,4]
[1,] 1.00 0.70 0.59 0.67
[2,] 0.70 1.00 0.72 0.73
[3,] 0.59 0.72 1.00 0.74
[4,] 0.67 0.73 0.74 1.00

and with 100,000 samples the calculated correlation is:
> Y = rmvnorm(100000,mean=rep(0,4),sigma=sig,method="svd")
> round(cor(Y),2)
     [,1] [,2] [,3] [,4]
[1,]  1.0  0.6  0.6  0.6
[2,]  0.6  1.0  0.7  0.7
[3,]  0.6  0.7  1.0  0.7
[4,]  0.6  0.7  0.7  1.0




A: I don't know R, so here are verbal directions.
Step 1. Let $U$ be the upper Cholesky factor of the correlation matrix.
For the correlation matrix you give, $U$ is
1  .6  .6       .6
0  .8  .425     .425
0   0  .677772  .235145
0   0   0       .635674

Step 2. Let $X$ be an $n \times 4$ matrix in which the first column is your given vector R,
with the three other columns filled with random independent standard normals.
Step 3. Subtract its mean from each column of $X$, then get the QR decomposition of the result.
If $R_{1,1}$ is negative then change the signs of all four values in row 1 of $U$.
Step 4. Let $Y = Q\,U \sqrt{n-1}$. If all you want is standard scores then you're done. Otherwise you can replace the first column of $Y$ by your given vector R, then for each of the three other columns, multiply by the desired standard deviation and add the desired mean.
