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
}
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
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. Chris Guure On Saturday, August 1, 2015 3:32 AM, David Winsemius <dwinsemius at 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
I am not sure how you are doing this but there is a package on CRAN which implements the Copas model (metasens). I am not sure whether that would help in your modelling. On 01/08/2015 02:36, 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 > > > } > > ______________________________________________ > R-help at 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. >-- Michael http://www.dewey.myzen.co.uk/home.html