On Jul 31, 2015, at 8:41 PM, Christopher Kelvin wrote:
> Thanks Dave.
>
> What I actually want is to obtain say 10, different sets of (n=50) data for
> every 10,000 iterations I run. You will realise that the current code
> produces one set of data (n=50). I want 10 different sets of 50 observations
> at one run. I hope this makes sense.
I would think that either `replicate(50, ...)` or `for(i in 1:50) {...}` would
suffice. Unless of course the phrase "want 10 different sets of 50 observations
at one run` means something different than it appears to request.
>
> Chris Guure
>
>
>
> On Saturday, August 1, 2015 3:32 AM, David Winsemius <dwinsem...@comcast.net>
> wrote:
>
>
> 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
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.
