####  Sample size calculations for the Bayesian Hui-Walter conditional independence model. 

### Splus 2000 code that generates the data sets and analyzes them using the WinBUGS 1.4.1 code below. 


library(embed)
rmultinom <- function(n, pr, long=F) {
  k <- length(pr)
    if(long) {
    y <- runif(n, 0, 1)
    p <- cumsum(pr)
    Seq <- 1:n
    x <- sapply(y, function(y, Seq, p) {Seq[y <= p][1]}, Seq=Seq, p=p)
  } else {
    x <- rep(NA,k)
    p <- pr/c(1,(1-cumsum(pr[1:(k-1)])))
    for (i in 1:(k-1)) {
      if (n==0) {
        x[i] <- 0
        if (i==k-1) x[k] <- 0
        next
      }
      y <- rbinom(1,n,p[i])
      x[i] <- y
      if (i==k-1) x[k] <- n-y
      n <- n-y
    }
  }
  return(x)
}


n1star <- ## the values of n1 go here
n2star <- ## the values of n2 go here

gen.data.mult <- function(n1=n1star, n2=n2star, Se2=0.80, Sp2=1)                          
   {
	p1 <- rep(0,4)
	p2 <- rep(0,4)
	Se1 <- rbeta(1, 42.6, 5.6)  
	Sp1 <- rbeta(1, 33.4, 33.4)
	pi1 <- rbeta(1, 42, 28.3)  
   	pi2 <-  rbeta(1, 28.3, 42) 
   	p1[1] <- pi1*Se1*Se2         + (1-pi1)*(1-Sp1)*(1-Sp2)
   	p1[2] <- pi1*Se1*(1-Se2)     + (1-pi1)*(1-Sp1)*Sp2
   	p1[3] <- pi1*(1-Se1)*Se2     + (1-pi1)*Sp1*(1-Sp2)
   	p1[4] <- 1 - (p1[1] + p1[2] + p1[3])
   	p2[1] <- pi2*Se1*Se2         + (1-pi2)*(1-Sp1)*(1-Sp2)
   	p2[2] <- pi2*Se1*(1-Se2)     + (1-pi2)*(1-Sp1)*Sp2
   	p2[3] <- pi2*(1-Se1)*Se2     + (1-pi2)*Sp1*(1-Sp2)
   	p2[4] <- 1 - (p2[1] + p2[2] + p2[3])
   	y1 <- rmultinom(n1, c(p1[1],p1[2],p1[3],p1[4]))
	y2 <- rmultinom(n2, c(p2[1],p2[2],p2[3],p2[4]))
	return(y1, y2, n1, n2, Sp2)
   }

write.multifile(func=gen.data.mult, towhere="c:/Program Files/WinBUGS14/SpS1/CondInd/r", n=1000, do.warning=F)

write.scripts(masterprog="C:/Program Files/WinBUGS14/SpS1/CondInd/model.txt",
              scriptroot="C:/Program Files/WinBUGS14/SpS1/CondInd/script",
              dataroot="C:/Program Files/WinBUGS14/SpS1/CondInd/r",
              saveroot="SpS1/CondInd/out",              monitor.list=c("eta2Se"),
              initfile="C:/Program Files/WinBUGS14/SpS1/CondInd/init.txt",gen.inits=F,
              nsims=100, nloops=10, update=5000, beg=1000)

#run.scripts1 _ function (nloops = 10, dos.script.root, 
#                        dos.location = "c:\\progra~1\\winbug~2\\winbug~1", 
#                        is.windows.2000 = T) 
#{
#    for (i in 1:nloops) {
#        cat("Running script #", i, " \n", sep = "")
#        system(paste(ifelse(is.windows.2000,"","start /w /m "),dos.location, " /par ", 
#            dos.script.root, as.character(i), ".txt", sep = ""), 
#            minimized = T)
#    }
#}

run.scripts1(nloops=10, dos.script.root="SpS1\\CondInd\\script",
             dos.location="C:\\Program Files\\WinBUGS14\\WinBUGS14.exe")


read.textdump("eta2Se", fromwhereroot="c:/Program Files/WinBUGS14/SpS1/CondInd/out", 
                        towhere=paste("c:/Program Files/WinBUGS14/SpS1/CondInd/results1", ".txt", sep = ""),
                        nloops=10, nsims=100)

p1Segp0 <- read.table("c:/Program Files/WinBUGS14/SpS1/CondInd/results1.txt", sep="", skip=1)[,2]








### WinBUGS code for the Bayesian Hui-Walter model


model;
{
y1[1:4] ~ dmulti(p1[1:4], n1)
y2[1:4] ~ dmulti(p2[1:4], n2)
p1[1] <- pi1*Se1*Se2         + (1-pi1)*(1-Sp1)*(1-Sp2)
p1[2] <- pi1*Se1*(1-Se2)     + (1-pi1)*(1-Sp1)*Sp2
p1[3] <- pi1*(1-Se1)*Se2     + (1-pi1)*Sp1*(1-Sp2)
p1[4] <- pi1*(1-Se1)*(1-Se2) + (1-pi1)*Sp1*Sp2
p2[1] <- pi2*Se1*Se2         + (1-pi2)*(1-Sp1)*(1-Sp2)
p2[2] <- pi2*Se1*(1-Se2)     + (1-pi2)*(1-Sp1)*Sp2
p2[3] <- pi2*(1-Se1)*Se2     + (1-pi2)*Sp1*(1-Sp2)
p2[4] <- pi2*(1-Se1)*(1-Se2) + (1-pi2)*Sp1*Sp2
pi1 ~ dbeta(125.8, 2.26)
pi2 ~ dbeta(2.26, 125.8)
Se1 ~ dbeta(375.3,375.3)
Sp1 ~ dbeta(316.5,36.1)
Se2 ~ dunif(0.5,1)
eta2Se <- step(Se2-0.90)
}