#### valdes.r ####

## Introduction of the function ##

# The 'valdes()' function generates 10000 posterior samples for log odds ratio and other model parameters based on a given simulated dataset


## Input Variables ##
# case.cmp: 2 by 2 validation table for cases
# cntr.cmp: 2 by 2 validation table for controls
# case.sum: 2 by 2 table for the main data in case group
# cntr.sum: 2 by 2 table for the main data in control group
# PRI: prior knowledge of the type of misclassification

## Output data ##
# res: a 10000 by 7 matrix containing 10000 samples from the posterior for model parameters
# count.rej: a 10000 by 7 matrix containing indictors for rejection of the proposed value for model parameters in MCMC process


valdes <- function(case.cmp, cntr.cmp, case.sum, cntr.sum, PRI) {

#initial values
burnin<-1000
n.iter<-10000+burnin
rho.r<-0.3 ## rho(r1,r2)
rho.s<-NA  ## rho(SN1,SN2)=rho(SP1,SP2)
rho1<-0 ## diff
rho2<-0.95 ## nearly non-diff.

# thetai=log(ri/(1-ri)), i=1,2
# pi=log(SNi/(1-SNi))
# qi=log(SPi/(1-SPi))

theta1.mu<-theta2.mu<- -1.94591
p1.mu<-p2.mu<-q1.mu<-q2.mu<-1.674952

theta1.sd<-theta2.sd<-0.9928112
p1.sd<-q1.sd<-p2.sd<-q2.sd<-0.6476974


theta1<-theta2<-theta1.mu
p1<-p2<-p1.mu
q1<-q2<-q1.mu

r1<-exp(theta1)/(1+exp(theta1))
r2<-exp(theta2)/(1+exp(theta2))
lor<-log(r2)-log(1-r2)+log(1-r1)-log(r1)

SN1<-exp(p1)/(1+exp(p1))
SN2<-exp(p2)/(1+exp(p2))

SP1<-exp(q1)/(1+exp(q1))
SP2<-exp(q2)/(1+exp(q2))    

case.mis <- cntr.mis <- case.tot <- cntr.tot <- matrix(NA,2,2)
res <- matrix(NA,n.iter,7)
colnames(res)<-c("r1","r2","SP1","SP2","SN1","SN2","logOR")
res[1,]<-c(r1,r2,SP1,SP2,SN1,SN2,lor)
count.rej<-matrix(NA,n.iter,6)
colnames(count.rej)<-c("r1","r2","SP1","SP2","SN1","SN2")


for (j in 2:(n.iter)) {

a1.case<-SP2*(1-r2)/(SP2*(1-r2)+(1-SN2)*r2)
a2.case<-SN2*r2/(SN2*r2+(1-SP2)*(1-r2))
a1.cntr<-SP1*(1-r1)/(SP1*(1-r1)+(1-SN1)*r1)
a2.cntr<-SN1*r1/(SN1*r1+(1-SP1)*(1-r1))

### update case.mis table (col sums = case.sum)

### first col (apparently exposed)
tmp1 <- rbinom(1,size=case.sum[1], prob=a2.case)
case.mis[,1] <- c(tmp1, case.sum[1]-tmp1)

### second col (apparently unexposed)
tmp2 <- rbinom(1,size=case.sum[2],prob=a1.case)
case.mis[,2] <- c(case.sum[2]-tmp2, tmp2)

### update cntr.mis table (col sums = cntr.sum)

### first col (apparently exposed)
tmp3 <- rbinom(1,size=cntr.sum[1],prob=a2.cntr)
cntr.mis[,1] <- c(tmp3 , cntr.sum[1]-tmp3)

### second col (apparently unexposed)
tmp4 <- rbinom(1,size=cntr.sum[2],prob=a1.cntr)
cntr.mis[,2] <- c(cntr.sum[2]-tmp4, tmp4)

### complete actual/apparent exposure tables
case.tot <- case.cmp+case.mis
cntr.tot <- cntr.cmp+cntr.mis

### update cntr prevalence

r1.tmp <- rbeta(1,sum(cntr.tot[1,])+1,sum(cntr.tot[2,])+1)

ratio1<-(dnorm(log(r1.tmp)-log(1-r1.tmp),mean=theta1.mu+rho.r*(log(r2)-log(1-r2)-theta1.mu),
         sd=theta1.sd*sqrt(1-rho.r^2)) /r1.tmp/(1-r1.tmp)) / 
        (dnorm(log(r1)-log(1-r1),mean=theta1.mu+rho.r*(log(r2)-log(1-r2)-theta1.mu),
         sd=theta1.sd*sqrt(1-rho.r^2)) /r1/(1-r1))

u<-runif(1)
if (u<=min(ratio1,1)) {
r1<-r1.tmp
count.rej[j,1]<-0
} else count.rej[j,1]<-1

### update case prevalence

r2.tmp <- rbeta(1,sum(case.tot[1,])+1,sum(case.tot[2,])+1)
ratio2<-(dnorm(log(r2.tmp)-log(1-r2.tmp),mean=theta1.mu+rho.r*(log(r1)-log(1-r1)-theta1.mu),
         sd=theta1.sd*sqrt(1-rho.r^2)) /r2.tmp/(1-r2.tmp)) / 
        (dnorm(log(r2)-log(1-r2),mean=theta1.mu+rho.r*(log(r1)-log(1-r1)-theta1.mu),
         sd=theta1.sd*sqrt(1-rho.r^2)) /r2/(1-r2))

u<-runif(1)
if (u<=min(ratio2,1)) {
r2<-r2.tmp
count.rej[j,2]<-0
} else count.rej[j,2]<-1


### update sensitivity and specificity

if (PRI=="NDF") {
  ### assuming nondifferentiality
  SN1.tmp<-SN2.tmp<-rbeta(1,(case.tot+cntr.tot)[1,1]+1,(case.tot+cntr.tot)[1,2]+1)
  
  ratio3<-(dnorm(log(SN1.tmp)-log(1-SN1.tmp),mean=p1.mu,sd=p1.sd)/SN1.tmp/(1-SN1.tmp)) /
          (dnorm(log(SN1)-log(1-SN1),mean=p1.mu,sd=p1.sd)/SN1/(1-SN1))
  u<-runif(1)
  if (u<=min(ratio3,1)) {
  SN1<-SN1.tmp
  SN2<-SN2.tmp
  count.rej[j,3]<-count.rej[j,4]<-0
  } else count.rej[j,3]<-count.rej[j,4]<-1  

  SP1.tmp<-SP2.tmp<-rbeta(1,(case.tot+cntr.tot)[2,2]+1,(case.tot+cntr.tot)[2,1]+1)
  ratio5<-(dnorm(log(SP1.tmp)-log(1-SP1.tmp),mean=q1.mu,sd=q1.sd)/SP1.tmp/(1-SP1.tmp)) /
        (dnorm(log(SP1)-log(1-SP1),mean=q1.mu,sd=q1.sd)/SP1/(1-SP1))

  u<-runif(1)
  if (u<=min(ratio5,1)) {
  SP1<-SP1.tmp
  SP2<-SP2.tmp
  count.rej[j,5]<-count.rej[j,6]<-0
  } else count.rej[j,5]<-count.rej[j,6]<-1
} #end PRI=NDF


if (PRI=="DF") {
  ### assuming differentiality, ind normal priors 

  rho.s<-rho1
  SN1.tmp<- rbeta(1, cntr.tot[1,1]+1, cntr.tot[1,2]+1 )
  ratio3<-(dnorm(log(SN1.tmp)-log(1-SN1.tmp),mean=p1.mu+rho.s*(log(SN2)-log(1-SN2)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN1.tmp/(1-SN1.tmp)) / 
         (dnorm(log(SN1)-log(1-SN1),mean=p1.mu+rho.s*(log(SN2)-log(1-SN2)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN1/(1-SN1))


  u<-runif(1)
  if (u<=min(ratio3,1)) {
  SN1<-SN1.tmp
  count.rej[j,3]<-0
  } else count.rej[j,3]<-1

  SN2.tmp <- rbeta(1, case.tot[1,1]+1, case.tot[1,2]+1 )
  ratio4<-(dnorm(log(SN2.tmp)-log(1-SN2.tmp),mean=p1.mu+rho.s*(log(SN1)-log(1-SN1)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN2.tmp/(1-SN2.tmp)) / 
         (dnorm(log(SN2)-log(1-SN2),mean=p1.mu+rho.s*(log(SN1)-log(1-SN1)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN2/(1-SN2))

  u<-runif(1)
  if (u<=min(ratio4,1)) {
  SN2<-SN2.tmp
  count.rej[j,4]<-0
  } else count.rej[j,4]<-1
  
  SP1.tmp <- rbeta(1, cntr.tot[2,2]+1, cntr.tot[2,1]+1 ) 
  ratio5<-(dnorm(log(SP1.tmp)-log(1-SP1.tmp),mean=q1.mu+rho.s*(log(SP2)-log(1-SP2)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP1.tmp/(1-SP1.tmp)) / 
         (dnorm(log(SP1)-log(1-SP1),mean=q1.mu+rho.s*(log(SP2)-log(1-SP2)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP1/(1-SP1))

  u<-runif(1)
  if (u<=min(ratio5,1)) {
  SP1<-SP1.tmp
  count.rej[j,5]<-0
  } else count.rej[j,5]<-1
  
  SP2.tmp <- rbeta(1, case.tot[2,2]+1, case.tot[2,1]+1 )
  ratio6<-(dnorm(log(SP2.tmp)-log(1-SP2.tmp),mean=q1.mu+rho.s*(log(SP1)-log(1-SP1)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP2.tmp/(1-SP2.tmp)) / 
         (dnorm(log(SP2)-log(1-SP2),mean=q1.mu+rho.s*(log(SP1)-log(1-SP1)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP2/(1-SP2))

  u<-runif(1)
  if (u<=min(ratio6,1)) {
  SP2<-SP2.tmp
  count.rej[j,6]<-0
  } else count.rej[j,6]<-1
} #end PRI=DF


if (PRI=="NNDF") {
  ### assuming differentiality, ind unif priors for the irreducible kernel
  rho.s<-rho2
  SN1.tmp<- rbeta(1, cntr.tot[1,1]+1, cntr.tot[1,2]+1 )
  ratio3<-(dnorm(log(SN1.tmp)-log(1-SN1.tmp),mean=p1.mu+rho.s*(log(SN2)-log(1-SN2)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN1.tmp/(1-SN1.tmp)) / 
         (dnorm(log(SN1)-log(1-SN1),mean=p1.mu+rho.s*(log(SN2)-log(1-SN2)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN1/(1-SN1))
  u<-runif(1)
  if (u<=min(ratio3,1)) {
  SN1<-SN1.tmp
  count.rej[j,3]<-0
  } else count.rej[j,3]<-1

  SN2.tmp <- rbeta(1, case.tot[1,1]+1, case.tot[1,2]+1 )
  ratio4<-(dnorm(log(SN2.tmp)-log(1-SN2.tmp),mean=p1.mu+rho.s*(log(SN1)-log(1-SN1)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN2.tmp/(1-SN2.tmp)) / 
         (dnorm(log(SN2)-log(1-SN2),mean=p1.mu+rho.s*(log(SN1)-log(1-SN1)-p1.mu),
         sd=p1.sd*sqrt(1-rho.s^2)) /SN2/(1-SN2))

  u<-runif(1)
  if (u<=min(ratio4,1)) {
  SN2<-SN2.tmp
  count.rej[j,4]<-0
  } else count.rej[j,4]<-1
  
  SP1.tmp <- rbeta(1, cntr.tot[2,2]+1, cntr.tot[2,1]+1 ) 
  ratio5<-(dnorm(log(SP1.tmp)-log(1-SP1.tmp),mean=q1.mu+rho.s*(log(SP2)-log(1-SP2)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP1.tmp/(1-SP1.tmp)) / 
         (dnorm(log(SP1)-log(1-SP1),mean=q1.mu+rho.s*(log(SP2)-log(1-SP2)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP1/(1-SP1))

  u<-runif(1)
  if (u<=min(ratio5,1)) {
  SP1<-SP1.tmp
  count.rej[j,5]<-0
  } else count.rej[j,5]<-1
  
  SP2.tmp <- rbeta(1, case.tot[2,2]+1, case.tot[2,1]+1 )
  ratio6<-(dnorm(log(SP2.tmp)-log(1-SP2.tmp),mean=q1.mu+rho.s*(log(SP1)-log(1-SP1)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP2.tmp/(1-SP2.tmp)) / 
         (dnorm(log(SP2)-log(1-SP2),mean=q1.mu+rho.s*(log(SP1)-log(1-SP1)-q1.mu),
         sd=q1.sd*sqrt(1-rho.s^2)) /SP2/(1-SP2))

  u<-runif(1)
  if (u<=min(ratio6,1)) {
  SP2<-SP2.tmp
  count.rej[j,6]<-0
  } else count.rej[j,6]<-1
}# end PRI=NNDF


lor<-log(r2)-log(1-r2)+log(1-r1)-log(r1)
res[j,] <- c(r1,r2,SP1,SP2,SN1,SN2,lor)

} #end if j

return(list(res[(burnin+1):n.iter,],count.rej[(burnin+1):n.iter,]))
}