[R] Help with automation of WinBUGS from R

[Anamika Chaudhuri]

Hi:
I have a question related to the R Code which calls BUGS. I have run the
model in WinBUGS and it runs fine giving me the expected results.  Below is
the automation code used when I had single outcome or univariate data for
Y’s. Now I want to use it for multiple outcomes. So the trial.data below
had one subscript earlier and now its two. There are 3 simulations for
testing.Not sure where to add the subscript in the code so that the same
process can be repeated for 2 outcomes instead of one.
library("R2WinBUGS")

> #Set working directory
> setwd("C://Tina/USB_Backup_042213/Testing")
> trial.data <- read.table("MVN_Test.csv", header=T)
> trial.data
  X.V11.V21.V31.V41.V51.V12.V22.V32.V42.V52
1                1,11,6,8,5,25,13,1,13,8,22
2                 2,9,1,7,9,25,13,1,18,9,12
3               3,6,4,9,10,26,15,1,14,11,12
>
> p1_true<- read.table("p1.csv",header=F)
> p1_true
          V1
1 0.15171877
2 0.07220343
3 0.13898233
4 0.13129042
5 0.43438334
>
> p2_true<-read.table("p2.csv",header=F)
> p2_true
          V1
1 0.20562945
2 0.04400286
3 0.25019831
4 0.18854418
5 0.22946289
> bugs.output <- list()

for(i in 1:3){
       Y <- as.integer(trial.data[i,])
       bugs.output[[i]] <- bugs(
       data=list(Y=as.matrix(Y), Nf=20, n=60,), # how to pass the other
parameters like mn which is a vector, prec which is a matrix, R a matrix
here?,
 inits=list(

list(gamma=c(0,0),T=structure(.Data=c(0.5,0,0.5,0),.Dim=c(2,2)))
               ),
       model.file="M-LN_model_trial.txt",
       parameters.to.save = c("p","rho","sigma2"),
n.chains=1, n.iter=12000, n.burnin=5000,
bugs.directory="C://Tina/USB_Backup_042213/winbugs14/WinBUGS14",
working.directory=NULL)}
bugs.output[[i]]$sims.list$p_trans<-t(bugs.output[[i]]$sims.list$p)
bugs.output[[i]]$sims.list$coverage<-matrix(data=0,nrow=20,ncol=1)
bugs.output[[i]]$sims.list$int_width<-matrix(data=0,nrow=20,ncol=1)
bugs.output[[i]]$sims.list$sqr_err<-matrix(data=0,nrow=20,ncol=1)
bugs.output[[i]]$sims.list$abs_bias<-matrix(data=0,nrow=20,ncol=1)

bugs.output[[i]]$sims.list$p_true<-p_true
for (j in 1:5) # change for no of sites
{
if
# define coverage
((bugs.output[[i]]$sims.list$p_true[j,]>quantile(bugs.output[[i]]$sims.list$p_trans[j,],c(.025)))
&
(bugs.output[[i]]$sims.list$p_true[j,]<quantile(bugs.output[[i]]$sims.list$p_trans[j,],c(.975))))

bugs.output[[i]]$sims.list$coverage[j,]<-1
else
bugs.output[[i]]$sims.list$coverage[j,]<-0
#define interval width
bugs.output[[i]]$sims.list$int_width[j,]<-(quantile(bugs.output[[i]]$sims.list$p_trans[j,],c(.975)))-(quantile(bugs.output[[i]]$sims.list$p_trans[j,],c(.025)))
# define sqr error loss
bugs.output[[i]]$sims.list$sqr_err[j,]<-(((bugs.output[[i]]$sims.list$p_true[j,])-mean(bugs.output[[i]]$sims.list$p_trans[j,]))^2)/0.17
# change with overall rate

# define absolute bias
bugs.output[[i]]$sims.list$abs_bias[j,]<-abs((bugs.output[[i]]$sims.list$p_true[j,])-mean(bugs.output[[i]]$sims.list$p_trans[j,]))
}
}

cov<-cbind(NULL) # initialize variable matrix cov before loop
interval_width<-cbind(NULL) # initialize variable matrix cov before loop
sqr_error_loss<-cbind(NULL) # initialize variable matrix square error loss
absolute_bias<-cbind(NULL) # initialize variable matrix absolute bias
for(i in 1:300){
 cov<-cbind(cov, bugs.output[[i]]$sims.list$coverage)
interval_width<-cbind(interval_width, bugs.output[[i]]$sims.list$int_width)
sqr_error_loss<-cbind(sqr_error_loss, bugs.output[[i]]$sims.list$sqr_err)
absolute_bias<-cbind(absolute_bias, bugs.output[[i]]$sims.list$abs_bias)
}
# Average Coverage
avg_cov<-apply(cov,1,mean)
cov_prob<-mean(avg_cov)
1-cov_prob
# Average Interval Width
avg_int_width<-apply(interval_width,1,mean)
avg_interval_width<-mean(avg_int_width)
avg_interval_width

# Average Square error loss
avg_err_loss<-apply(sqr_error_loss,1,mean)
average_loss<-mean(avg_err_loss)
average_loss

# Average abosulte Bias
avg_abs_bias<-apply(absolute_bias,1,mean)
average_bias<-mean(avg_abs_bias)
average_bias


Thanks for your help.

Anamika

        [[alternative HTML version deleted]]


______________________________________________
R-help@r-project.org mailing list
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.