pba.tri<-function(a,b,c,d,
                  se1.a,se1.b,se1.c,
                  sp1.a,sp1.b,sp1.c,
                  se0.a,se0.b,se0.c,
                  sp0.a,sp0.b,sp0.c,
                  nd,SIMS){
  
  n_case <- a + b
  n_ctrl <- c + d
  
  # set some important parameters
  niter <- SIMS #number of iterations
  
  # draw sensitivities and specificities
  if(nd==0){
    se1<-rtriangle(niter,se1.a,se1.c,se1.b)
    sp1<-rtriangle(niter,sp1.a,sp1.c,sp1.b)
    se0<-rtriangle(niter,se0.a,se0.c,se0.b)
    sp0<-rtriangle(niter,sp0.a,sp0.c,sp0.b)
  }
  if(nd==1){
    se0<-se1<-rtriangle(niter,se1.a,se1.c,se1.b)
    sp0<-sp1<-rtriangle(niter,sp1.a,sp1.c,sp1.b)
  }
  
  
  # calculate bias-adjusted cell frequencies
  ac <- (a - n_case*(1-sp1))/(se1 - (1-sp1)) #bias-adjusted cases, exposed 
  bc <- n_case - ac #bias-adjusted cases, unexposed
  cc <- (c - n_ctrl*(1-sp0))/(se0 - (1-sp0)) #bias-adjusted controls, exposed
  dc <- n_ctrl - cc #bias-adjusted controls, unexposed
  flag<-(ac<=0)|(bc<=0)|(cc<=0)|(dc<=0)
  
  #calculate bias-adjusted odds ratio, second source of uncertainty
  #calculate prevalence of exposure in cases and controls, accounting for sampling error
  PrevE_cases<-rbeta(niter,ac,bc)
  PrevE_controls <- rbeta(niter,cc,dc)
  
  #calculate PPV and NPV of exposure classification in cases and controls
  PPV_case <- (se1*PrevE_cases)/((se1*PrevE_cases)+(1-sp1)*(1-PrevE_cases))
  PPV_control <- (se0*PrevE_controls)/((se0*PrevE_controls)+(1-sp0)*(1-PrevE_controls))
  NPV_case <- (sp1*(1-PrevE_cases))/((1-se1)*PrevE_cases+sp1*(1-PrevE_cases))
  NPV_control <- (sp0*(1-PrevE_controls))/((1-se0)*PrevE_controls+sp0*(1-PrevE_controls))
  
  #calculate the expected number of exposed cases and controls
  ab <- rbinom(niter,a,PPV_case)+rbinom(niter,b,1-NPV_case)
  bb <- n_case-ab
  cb <- rbinom(niter,c,PPV_control)+rbinom(niter,d,1-NPV_control)
  db <- n_ctrl-cb
  
  #calculate bias adjusted OR with second source of uncertainty
  or_bb <- (ab/cb)/(bb/db)
  
  # calculate bias-adjusted odds ratios, third source of uncertainty, bias-adjusted standard error
  se_bb <- sqrt(1/ab+1/bb+1/cb+1/db)
  z <- rnorm(niter)
  or_bb_cb <- exp(log(or_bb) - (z*se_bb))
  
  #or resampel prev to introduce random noise
  pi1=rbeta(niter,ab,bb)
  pi0=rbeta(niter,cb,db)
  or_1=(pi1/(1-pi1))/(pi0/(1-pi0))
  
  # bind vectors as dataframe
  bound <- data.frame(ac, bc, cc, dc, ab, bb, cb, db, or_bb_cb, or_1)
  
  # housekeeping - retain rows with positive corrected cell freqs
  summary_pba <- bound %>%
    filter(ac > 0) %>%
    filter(bc > 0) %>%
    filter(cc > 0) %>%
    filter(dc > 0) %>%
    filter(ab > 0) %>%
    filter(bb > 0) %>%
    filter(cb > 0) %>%
    filter(db > 0) 
  
  q3<-quantile(summary_pba$or_bb_cb,c(0.025,0.5,0.975))
  q4<-quantile(summary_pba$or_1,c(.025,.5,.975))
  
  return(list(c("newABCD"=q3), c("OR"=summary_pba$or_bb_cb)))
  
}