Hello, I am trying to do a simulation process, which consists of making a function to generate data, then fit a correlated model and another uncorrelated model with the generated data, this must be repeated a large number of times, for example 100. With the The code shown below only generates one output and I want it to generate the outputs or summaries of all the models adjusted simultaneously, that is, if there are 100 simulations, they should generate 200 summaries, since they are two models. Can you help me, please.

simul=function(var.v = 1.7, var.w = 2.6, corr = 0, n = 217, mean=52.5576, sd= 20.45235, mu=c(0,0), a0 = c(1.9, -0.01), b0 = c(2.5, -0.02)){

cov = sqrt(var.v*var.w) corr
sigma = matrix(c(var.v, cov, cov, var.w), 2, 2)
library(MASS)
u = mvrnorm(n = n, mu, Sigma=sigma)
x = rnorm(n, mean, sd)
X = cbind( rep(1, length(x)), x)
mu1 = exp(X%%a0 + u[,1])
mu2 = exp(X%*%b0 + u[,2])

y1 = rpois(n, mu1)

y2 = rpois(n, mu2)

y = c(y1, y2)

y

Xb = c(x, x)

id = c(1:n, 1:n)

time = as.factor(c(rep(1, n), rep(2, n)))

# plot(y1, y2)

datos=data.frame(y,Xb,time,id)

library(MASS)

library(lme4)

fit1=glm(y ~ Xb*time)
fit2<- glmer(y~Xb*time + (time-1 | id) , family = poisson, data=datos)

summary(fit1)

summary(fit2)

}

simul(var.v = 1.7, var.w = 2.6, corr = 0, n = 217, mean=52.5576, sd= 20.45235, mu=c(0,0), a0 = c(1.9, -0.01), b0 = c(2.5, -0.02))

sims = replicate(100, simul(), simplify = FALSE )

sims

**strong text**