On Jul 31, 2015, at 6:36 PM, Christopher Kelvin via R-help wrote:

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

OK. The code runs without error. So .... what exactly is the problem? (I have 
no experience with the Copas selection model if in fact that is what is being 
exemplified.)

-- 

David Winsemius
Alameda, CA, USA

______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to