[R] Clarification on Simulation and Iteration

[Christopher Kelvin]

Dear All,
I am performing some simulations for a new model. I run about 10,000 iterations with a sample of 50 datasets and this returns one set of 50 simulated data. 

Now, what I need to obtain is 10 sets of the 50 simulated data out of the 10,000 iterations and not just only 1 set.  The model is the Copas selection model for publication bias in Mete-analysis. Any one who knows this model has any suggestion for the improvement of my code is most welcome.

Below is my code. 


Kind regards


Chris Guure
University of Ghana




install.packages("msm") 
library(msm) 


rho1=-0.3; tua=0.020; n=50; d=-0.2; rr=10000; a=-1.3; b=0.06 
si<-rtnorm(n, mean=0, sd=1, lower=0, upper=0.2)# I used this to generate standard errors for each study 
set.seed(21111)   ## I have stored the data and the output in this seed 

for( i in 1:rr){ 

mu<-rnorm(n,d,tua^2)              # prob. of each effect estimate 
rho<-si*rho1/sqrt(tua^2 + si^2) # estimate of the correlation coefficient 
mu0<- a + b/si       # mean of the truncated normal model (Copas selection model) 
y1<-rnorm(mu,si^2)            # observed effects zise 
z<-rnorm(mu0,1)               # selection model 
rho2<-cor(y1, z) 

select<-pnorm((mu0 + rho*(y1-mu)/sqrt(tua^2 + si^2))/sqrt(1-rho^2)) 
probselect<-ifelse(select<z, y1, NA)# the prob that the study is selected 

probselect 
data<-data.frame(probselect,si)    # this contains both include and exclude data 
data 
data1<-data[complete.cases(data),] # Contains only the included data for analysis 
data1 


}