How to simulate data conditional on variables and respecting correlation structure in R [closed]

I am trying to simulate data to benchmark a multi-dimensional data method (referred to as a multi-omics approach) with a specific correlation structure and multiple confounding variables. For the first aspect I used a well-known algorithm to transform my data distributions called NORTA (Normal to Anything). Intuitively the NORTA algorithm allows one to specify the correlation between normal random variables and convert, through a quantile-quantile transformation, into any arbitrary statistical distribution. In my case normal to a zero-inflated negative binomial distribution.

I simulated 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

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. Is it possible to use the NORTA algorithm (or another statistical transformation) to simulate non-normal data with a specific correlation structure ? The NORTA algorithm seems to work well for standard multivariate normal distributions, but not for normal distribution with non-zero means. Indeed, if I first generate 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.

• Close-voters: I believe this question is primarily about statistics, not programming, and that it should therefore stay open/be reopened. See Closing "software questions". Mar 28, 2023 at 6:45

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:

1. Simulate a large multivariate normal data set with the correlations you want. I use N=1,000,000.
2. Transform each variable to your desired distribution.
3. Compute the new correlation matrix, finding attenuated correlations.
4. For each correlation, do a regression between the point [0,0] and the point [attenuated correlation, simulated correlation].
5. Solve for what simulated correlation will result in your desired outcome correlation.
6. Update the correlations you're using to generate your data.
7. Compute some error (i.e., the difference between outcome and desired).
8. Simulate a new large correlation matrix, now using your updated correlations.
9. Repeat 2-7 until your error is within some tolerance, or until you reach some maximum number of iterations.
10. 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
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
sim.settings[ind]=temp.dat$$sim.settings ind2=ind ind2[,1]=ind[,2] ind2[,2]=ind[,1] sim.settings[ind2]=temp.dat$$sim.settings