# Generate a random variable with a defined correlation to multiple existing variables

this question is strongly related to: Generate a random variable with a defined correlation to an existing variable.

However I'm struggling to implement it in a more complex matter:

Given X, a covariance matrix of two variables f1, f2. I simulate data for f1 and f2, with the mvrnorm function from the MASS package. Now I want to simulate data for y1, y2, y3, y4 on the predefined correlation matrix with f1, f2 called Y. So to be more precise, yi has to correlate on multiple existing variables, in this case f1 and f2. So I wrote some code based on the first answer on this question Generate a random variable with a defined correlation to an existing variable. But as you can see, in the result matrix Y, I don't have an idea of how to simulate a correlation on multiple existing variables. Is it even possible with this method? I strongly appreciate hints, suggestions or readings.

Here's the code:

f1 <- c(1, 0.5)
f2 <- c(0.5, 1)

sig <- rbind(f1, f2)

colnames(sig) <- rownames(sig)

fscores <- mvrnorm(100, Sigma = sig, mu = c(0, 0), empirical = TRUE)

X <- cor(fscores)

fnames <- c(rep("f1", 2), rep("f2", 2))

fycoeffs <- c(0.7, 0.7, 0.7, 0.7)

sim.data <- matrix(numeric(), 100, length(fycoeffs))

for(i in seq_along(fycoeffs)){

sim.data[, i] <- sim_data(fycoeffs[i], fscores[, fnames[i]])
}

Y <- cor(sim.data, fscores)

rownames(Y) <- c("y1", "y2", "y3", "y4")

print(Y)

f1        f2
y1 0.7000000 0.3337438
y2 0.7000000 0.3309644
y3 0.7000000 0.3457498
y4 0.2754015 0.7000000
y5 0.2665275 0.7000000
y6 0.2959004 0.7000000

#' Simulate random variable with correlation rho on existing variable fscore.
#'
sim_data <- function(rho, fscore){

if(rho == 1){

result <- fscore
}else{

theta <- acos(rho)             # corresponding angle
x2    <- rnorm(100, 0, 1)      # new random data
X     <- cbind(fscore, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(100)                               # 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 fscore
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to Xctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

result <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
}

return(result)
}