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

Question

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

Answer 1

I think your question is more sort of philosophical than mathematical. Once you have a fixed number, it's a constant and does not correlate to anything, by the definition. But of course you can cheat and view the result as 'random' so that you just think that you get to know f1 and f2 before, and you want to merely simulate from the conditional distribution of yi given f. (Think about it. This is an important distinction.)

You can create correlated random variables in many ways. Wikipedia provides the conditional distribution of multivariate normal variables assuming that part of the components are known. I.e. you view all variables as one vector x = (f;y), and ask what's the conditional distribution of yi given f. If x is initially MVN, also (yi|f) will be.