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@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
______________________________________________
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.
