How to simulate data conditional on variables and respecting correlation structure in R I try to simulate data for a benchmark multi-dimensional data methods (refer to as multi-omics approaches) with specific correlation structures and depending on others variables. For the first aspect I use a well-known algorithm called NORTA (Normal to Anything). Intuitively the NORTA algorithm allows to specify correlation between normal random variables and convert through quantile-quantile transformation into any arbitrary statistical distributions. In my case normal to zero-inflated negative binomial distribution.
In my case, I simulate 25 correlated variables using NORTA algorithm like this:
set.seed(1234)
L0 = matrix(0, ncol=25, nrow=25)
#Variance for microbiotes throughout samples
#We can play on this parameter in order to assess biological variability
diag(L0) = runif(25,1.5,2.5)
  
#Off-diagonal elements are randomly selected to have either 0 covariance or a positive or negative covariance based on uniform distribution
L0[lower.tri(L0)] = sapply(1:length(L0[lower.tri(L0)]), function(x) sample(c(0,runif(1,-1.5,1.5)),1, prob = c(0.7,0.3)))
  
Precision0 = L0%*%t(L0)
  
#We obtain Covariance matrix based on Cholesky decomposition of lower triangular matrix
Sigma0 = solve(Precision0)
Cor0 = cov2cor(Sigma0)
  
#The multivariate normal distribution is generated for 100 individuals with mean 0 and the 
#Correlation structure 
multi.norm = MASS::mvrnorm(100, rep(0,25), Cor0)  

#Now we can simulate data from zero-inflated distribution while conserving original correlation structure
simulated.microbiotes = matrix(VGAM::qzinegbin(pnorm(multi.norm), size=0.3763196,mu=exp(10.12693), pstr0 = 0.3),ncol=25, nrow=100)

Then I am interested in simulating metabolite data depending on certain microbiotes. For one microbiote impacting one metabolite I think I can generate associated data, assuming a mean depending on my microbiote level multiply by a certain coefficient.
random.microbiote = sample(1:ncol(microbiotes), 1)
coef = 0.8

#Here I used the MVNORM object since using the Zero-inflated negative binomial distribution
#leads to incorrect results

simulated.metabolite = VGAM::rzinegbin(100, munb=exp(10.12693 + coef *multi.norm[,random.microbiote]), size=0.38)

Using replicate function I can easily simulate 50 variables where X% depend on Y% of my microbiotes, assuming different association levels. However, my variables will be correlated only by chance. I would like to specify certain levels of correlation, by hand. I wonder if I can use the NORTA algorithm to simulate data depending on microbiotes while assuming correlation structure ? NORTA algorithm seems to work well for standard multivariate normal distribution not for normal distribution with non-zero means. Indeed, if I firstly generated correlated data and then change certain variables by changing their means depending on one specific microbiote, I destroy the original correlation structure, which is not wanted in my case.
I don't know if I can do this kind of thing easily, but any insights will be welcome.

Yes, NORTA attenuates correlations. You can use it to transform multiple variables so that they have the right distribution and correlation, but the tradeoff is you lose control of the multivariate symmetry. Depending on your goals, this may or may not be worth it.
The trick is to simulate a higher correlation so that the transformation results in the correct correlation. The general algorithm is as follows:

*

*Simulate a large multivariate normal data set with the correlations you want. I use N=1,000,000.

*Transform each variable to your desired distribution.

*Compute the new correlation matrix, finding attenuated correlations.

*For each correlation, do a regression between the point [0,0] and the point [attenuated correlation, simulated correlation].

*Solve for what simulated correlation will result in your desired outcome correlation.

*Update the correlations you're using to generate your data.

*Compute some error (i.e., the difference between outcome and desired).

*Simulate a new large correlation matrix, now using your updated correlations.

*Repeat 2-7 until your error is within some tolerance, or until you reach some maximum number of iterations.

*Use the final set of correlations to simulate a new data set you'll use for your analyses.

Edit: Here is example code using a binomial transformation, as it is a lot faster. It converges in 2 iterations. Note that you have to simulate a big data set, N=1 million, for this to work. It's effectively a brute-force strategy.
set.seed(1234)

var.num=10
obs.num=1000000 # Needs to be large to work - brute force strategy
tol = 0.005 #repeat until error is within this tolerance. 



L0 = matrix(0, ncol=var.num, nrow=var.num)
diag(L0) = runif(var.num,1.5,2.5)
L0[lower.tri(L0)] = sapply(1:length(L0[lower.tri(L0)]), function(x) sample(c(0,runif(1,-1.5,1.5)),1, prob = c(0.7,0.3)))
Precision0 = L0%*%t(L0)
Sigma0 = solve(Precision0)
Cor0 = cov2cor(Sigma0)

# this is the correlation matrix we are trying to approximate
target.cor = Cor0
# start simulations with target matrix
sim.settings = Cor0

# initalize error
error = 1

# while the error is greater than tol
while(error > tol){

  #simulate data
  multi.norm = MASS::mvrnorm(obs.num, rep(0,var.num), sim.settings)  
  # use a faster transformation for demonstration purposes
  simulated.microbiotes = matrix(stats::qbinom(pnorm(multi.norm), size = 200,prob = .01),ncol=var.num, nrow=obs.num)
  
  # correlation after transformation
  new.cor = cor(simulated.microbiotes)
  
  # compute error
  error = max(abs(target.cor[lower.tri(target.cor)] - new.cor[lower.tri(new.cor)]))
  print(error)
  
    temp.dat=data.frame(target = target.cor[lower.tri(target.cor)],
                       sim.settings = sim.settings[lower.tri(sim.settings)],
                       new.cor = new.cor[lower.tri(new.cor)])
  
  # estimate what simulation settings will yield the target after transformation
    adjust<-apply(temp.dat,MARGIN = 1,FUN = function(x){
      d<-base::data.frame(y =c(0,x[2]),x= c(0,x[3]))
      f1<-stats::lm( y~x,data = d )
      p0<-polynom::polynomial(stats::coefficients(f1))
      return(stats::predict(p0,x[1]))
    })
  
   # update simulation settings
    temp.dat$sim.settings = adjust
    
    # re-build the simulation setting matrix
    ind <- base::which( base::lower.tri(sim.settings,diag=F) , arr.ind = TRUE )
    
    sim.settings[ind]=temp.dat$sim.settings
    ind2=ind
    ind2[,1]=ind[,2]
    ind2[,2]=ind[,1]
    sim.settings[ind2]=temp.dat$sim.settings
    
}
 

