##### BayesMisc.r #####

## Introduction of the function ##

# The 'BayesMisc()' function calls another function 'valdes()' to estimate the log odds ratio (and other parameters) on simulated datasets over different scenarios 
# given the type of misclassification in prior belief.
# It produces the posterior mean and 100(1-alpha)% credible interval for log odds ratio for each dataset in a scenario.
# It creates mean-squared error of point estimators for log odds ratio, coverage and average width of nominal 95% interval estimators.



## Input Variables ##

# NSCEN: number of misclassification scenarios in simulation study (NSCEN=4 by default)
# NREP: number of simulated datasets in each scenario that will be analyzed using the Bayesian method (NREP=400 by default)
# alpha: type I error
# misc: degree of differentiality of misclassification in prior belief (misc="All" by default)
#       "DF" - differential misclassification
#       "NDF" - nondifferential misclassification
#       "NNDF" - nearly nondifferential misclassification
#       "All" - results will be generated for all differentiality levels above 
# my.tab: a preloaded data structure containing all simulated datasets


## Output Files ##

# Output data are saved in text files for each misclassification level ("DF", "NDF", "NNDF" and "All").
# A series of text files are generated for each misclassification level.

#1. "summary.txt":
#This file contains mean-squared error of point estimators for log odds ratio, coverage and average width of nominal 95% interval 
#estimators over NREP(400 by default) simulated datasets for each scenario

#2. "result1.txt", "result2.txt", "result3.txt", "result4.txt":
#One file is generated for each scenario.
#Each file contains the posterior mean and 100(1-alpha)% credible interval for logOR estimator, calculated from a sinle simulated 
#dataset in a given scenario (number of rows=NREP in each file).

#3. "avepar1.txt", "avepar2.txt", "avepar3.txt", "avepar4.txt":
#One file is generated for each scenario.
#Each file contains the posterior point estimator vector for model paramenters "r1","r2","sp1","sp2","sn1","sn2","lor" 
#(2 indicates case, 1 indicates control), calculated from a sinle simulated dataset in a given scenario (number of rows=NREP in each file).

#4. "averejec1.txt", "averejec2.txt", "averejec3.txt", "averejec4.txt":
#One file is generated for each scenario.
#Each file contains the percent of proposal rejection during sampling from the posterior for model parameters 
#"r1","r2","sp1","sp2","sn1","sn2","lor" (2 indicates case, 1 indicates control), calculated from a sinle simulated dataset 
#in a given scenario (number of rows=NREP in each file).
        
#5. "lastpar1.txt", "lastpar2.txt","lastpar3.txt","lastpar4.txt":
#One file is generated for each scenario.
#Each file contains 10000 posterior samples of model parameters "r1","r2","sp1","sp2","sn1","sn2","lor" based on 
#the last (the NREPth) simulated dataset in a given scenario.
   


BayesMisc<- function(NSCEN=4,NREP=400,alpha=0.05,misc="All") {

source("valdes.R")

#############################################
if (misc=="DF") {
# to invoke differential misclassification procedure
print("misc==DF")
summary <- matrix(NA,nrow=3*NSCEN,ncol=1)
colnames(summary)<-"DF"
rownames(summary)<-rep(c("MSE","Coverage Proportion","Ave. Width"),NSCEN)


for (SCEN in 1:NSCEN) {

cat("SCEN:",SCEN,"\n")

result <- matrix(NA,nrow=NREP,ncol=3)
colnames(result)<-rep(c("mean","lower bound","upper bound"),1)

sumr <- matrix(NA,nrow=3,ncol=1)
ave.param<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection)<-c("r1","r2","sp1","sp2","sn1","sn2")

for (i in 1:NREP) {

cat("i:",i,"\n")

table.tmp<-my.tab[,(2*i-1):(2*i),SCEN]
cntr.cmp<-table.tmp[1:2,]
case.cmp<-table.tmp[3:4,]
cntr.sum<-table.tmp[5,]
case.sum<-table.tmp[6,]

# Call function 'valdes()' to sample model parameters from the posterior distribution
rest.valdes.d<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"DF")
param<-rest.valdes.d[[1]]
ave.param[i,]<-apply(param,2,mean)

# rejection vectors: 0=accept proposal, 1=reject proposal
rejection<-rest.valdes.d[[2]]
ave.rejection[i,]<-apply(rejection,2,mean)

# mean and 100(1-alpha)% interval of logOR
lor<-param[,7] # extract logOR
result[i,] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))

} #end for i


flag1 <- ifelse((Logodd>=result[,2] & Logodd<=result[,3]),1,0)


# MSE
sumr[1,]<-mean((result[,1]-Logodd)^2)
# Coverage proportion
sumr[2,]<-sum(flag1)/NREP
# Average length of the 100*(1-alpha)% interval
sumr[3,]<-mean(result[,3]-result[,2])

summary[(3*(SCEN-1)+1):(3*SCEN),]<- sumr
print(result)

if (SCEN==1) {
        write.table(result,file="Bayes/DF/result1.txt",sep=" ")
        write.table(ave.param,file="Bayes/DF/avepar1.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/DF/averejec1.txt",sep=" ")
        # save param for the NREPth dataset in each scenario (i.e. the 400th dataset in each scenario, by default)
    write.table(param,file="Bayes/DF/lastpar1.txt",sep=" ")
}

if (SCEN==2) {
        write.table(result,file="Bayes/DF/result2.txt",sep=" ")
        write.table(ave.param,file="Bayes/DF/avepar2.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/DF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/DF/lastpar2.txt",sep=" ")
}

if (SCEN==3) {
        write.table(result,file="Bayes/DF/result3.txt",sep=" ")
        write.table(ave.param,file="Bayes/DF/avepar3.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/DF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/DF/lastpar3.txt",sep=" ")
}

if (SCEN==4) {
        write.table(result,file="Bayes/DF/result4.txt",sep=" ")
        write.table(ave.param,file="Bayes/DF/avepar4.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/DF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/DF/lastpar4.txt",sep=" ")
}


} #end for SCEN

print(summary)

write.table(summary,file="Bayes/DF/summary.txt",sep=" ")

} #end misc=DF

######################################################
if (misc=="NDF") {
# to invoke non-differential misclassification
print("misc==NDF")
summary <- matrix(NA,nrow=3*NSCEN,ncol=1)
colnames(summary)<-"NDF"
rownames(summary)<-rep(c("MSE","Coverage Proportion","Ave. Width"),NSCEN)


for (SCEN in 1:NSCEN) {

cat("SCEN:",SCEN,"\n")

result <- matrix(NA,nrow=NREP,ncol=3)
colnames(result)<-rep(c("mean","lower bound","upper bound"),1)

sumr <- matrix(NA,nrow=3,ncol=1)
ave.param<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection)<-c("r1","r2","sp1","sp2","sn1","sn2")

for (i in 1:NREP) {

cat("i:",i,"\n")

table.tmp<-my.tab[,(2*i-1):(2*i),SCEN]
cntr.cmp<-table.tmp[1:2,]
case.cmp<-table.tmp[3:4,]
cntr.sum<-table.tmp[5,]
case.sum<-table.tmp[6,]

rest.valdes.ndf<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"NDF")
param<-rest.valdes.ndf[[1]] 
ave.param[i,]<-apply(param,2,mean)

# rejection vectors: 0=accept, 1=reject
rejection<-rest.valdes.ndf[[2]]
ave.rejection[i,]<-apply(rejection,2,mean)

# mean and 100(1-alpha)% interval of logOR
lor<-param[,7] # extract logOR
result[i,] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))

} #end for i


flag1 <- ifelse((Logodd>=result[,2] & Logodd<=result[,3]),1,0)


# MSE
sumr[1,]<-mean((result[,1]-Logodd)^2)
# Coverage proportion
sumr[2,]<-sum(flag1)/NREP
# Average length of the 100(1-alpha)% interval
sumr[3,]<-mean(result[,3]-result[,2])

summary[(3*(SCEN-1)+1):(3*SCEN),]<- sumr
print(result)

if (SCEN==1) {
        write.table(result,file="Bayes/NDF/result1.txt",sep=" ")
        write.table(ave.param,file="Bayes/NDF/avepar1.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NDF/averejec1.txt",sep=" ")
        # save param for the NREPth dataset in each scenario (i.e. the 400th dataset in each scenario, by default)
    write.table(param,file="Bayes/NDF/lastpar1.txt",sep=" ")
}

if (SCEN==2) {
        write.table(result,file="Bayes/NDF/result2.txt",sep=" ")
        write.table(ave.param,file="Bayes/NDF/avepar2.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NDF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NDF/lastpar2.txt",sep=" ")
}

if (SCEN==3) {
        write.table(result,file="Bayes/NDF/result3.txt",sep=" ")
        write.table(ave.param,file="Bayes/NDF/avepar3.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NDF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NDF/lastpar3.txt",sep=" ")
}

if (SCEN==4) {
        write.table(result,file="Bayes/NDF/result4.txt",sep=" ")
        write.table(ave.param,file="Bayes/NDF/avepar4.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NDF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NDF/lastpar4.txt",sep=" ")
}


} #end for SCEN

print(summary)

write.table(summary,file="Bayes/NDF/summary.txt",sep=" ")

} #end misc=NDF


##########################################################
if (misc=="NNDF") {
# to invoke nearly non-differential misclassification
print("misc==NNDF")
summary <- matrix(NA,nrow=3*NSCEN,ncol=1)
colnames(summary)<-"NNDF"
rownames(summary)<-rep(c("MSE","Coverage Proportion","Ave. Width"),NSCEN)


for (SCEN in 1:NSCEN) {

cat("SCEN:",SCEN,"\n")

result <- matrix(NA,nrow=NREP,ncol=3)
colnames(result)<-rep(c("mean","lower bound","upper bound"),1)

sumr <- matrix(NA,nrow=3,ncol=1)
ave.param<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection)<-c("r1","r2","sp1","sp2","sn1","sn2")

for (i in 1:NREP) {

cat("i:",i,"\n")

table.tmp<-my.tab[,(2*i-1):(2*i),SCEN]
cntr.cmp<-table.tmp[1:2,]
case.cmp<-table.tmp[3:4,]
cntr.sum<-table.tmp[5,]
case.sum<-table.tmp[6,]

rest.valdes.nndf<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"NNDF")
param<-rest.valdes.nndf[[1]] 
ave.param[i,]<-apply(param,2,mean)

# rejection vectors: 0=accept, 1=reject
rejection<-rest.valdes.nndf[[2]]
ave.rejection[i,]<-apply(rejection,2,mean)

# mean and 100(1-alpha)% interval of logOR
lor<-param[,7] # extract logOR
result[i,] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))

} #end for i


flag1 <- ifelse((Logodd>=result[,2] & Logodd<=result[,3]),1,0)


# MSE
sumr[1,]<-mean((result[,1]-Logodd)^2)
# Coverage proportion
sumr[2,]<-sum(flag1)/NREP
# Average length of the 90% interval
sumr[3,]<-mean(result[,3]-result[,2])

summary[(3*(SCEN-1)+1):(3*SCEN),]<- sumr
print(result)

if (SCEN==1) {
        write.table(result,file="Bayes/NNDF/result1.txt",sep=" ")
        write.table(ave.param,file="Bayes/NNDF/avepar1.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NNDF/averejec1.txt",sep=" ")
       # save param for the NREPth dataset in each scenario (i.e. the 400th dataset in each scenario, by default)
    write.table(param,file="Bayes/NNDF/lastpar1.txt",sep=" ")
}

if (SCEN==2) {
        write.table(result,file="Bayes/NNDF/result2.txt",sep=" ")
        write.table(ave.param,file="Bayes/NNDF/avepar2.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NNDF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NNDF/lastpar2.txt",sep=" ")
}

if (SCEN==3) {
        write.table(result,file="Bayes/NNDF/result3.txt",sep=" ")
        write.table(ave.param,file="Bayes/NNDF/avepar3.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NNDF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NNDF/lastpar3.txt",sep=" ")
}

if (SCEN==4) {
        write.table(result,file="Bayes/NNDF/result4.txt",sep=" ")
        write.table(ave.param,file="Bayes/NNDF/avepar4.txt",sep=" ")
        write.table(ave.rejection,file="Bayes/NNDF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param,file="Bayes/NNDF/lastpar4.txt",sep=" ")
}


} #end for SCEN

print(summary)

write.table(summary,file="Bayes/NNDF/summary.txt",sep=" ")

} #end misc=NNDF


#######################################################

if (misc=="All") {
# to invoke all three routines 

print("misc==All")

summary <- matrix(NA,nrow=3*NSCEN,ncol=3)
colnames(summary)<-c("NDF","NNDF","DF")
rownames(summary)<-rep(c("MSE","Coverage Proportion","Ave. Width"),NSCEN)

for (SCEN in 1:NSCEN) {

cat("SCEN:",SCEN,"\n")

result <- matrix(NA,nrow=NREP,ncol=9)
colnames(result)<-rep(c("mean","lower bound","upper bound"),3)

sumr <- matrix(NA,nrow=3,ncol=3)

ave.param.ndf<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param.ndf)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection.ndf<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection.ndf)<-c("r1","r2","sp1","sp2","sn1","sn2")

ave.param.nndf<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param.nndf)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection.nndf<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection.nndf)<-c("r1","r2","sp1","sp2","sn1","sn2")

ave.param.df<-matrix(NA,nrow=NREP,ncol=7)
colnames(ave.param.df)<-c("r1","r2","sp1","sp2","sn1","sn2","lor")
ave.rejection.df<-matrix(NA,nrow=NREP,ncol=6)
colnames(ave.rejection.df)<-c("r1","r2","sp1","sp2","sn1","sn2")


for (i in 1:NREP) {

cat("i:",i,"\n")

table.tmp<-my.tab[,(2*i-1):(2*i),SCEN]
cntr.cmp<-table.tmp[1:2,]
case.cmp<-table.tmp[3:4,]
cntr.sum<-table.tmp[5,]
case.sum<-table.tmp[6,]

rest.valdes.ndf<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"NDF")
rest.valdes.nndf<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"NNDF")
rest.valdes.df<- valdes(case.cmp, cntr.cmp, case.sum, cntr.sum,"DF")


param.ndf<-rest.valdes.ndf[[1]] 
ave.param.ndf[i,]<-apply(param.ndf,2,mean)
# rejection vectors: 0=accept, 1=reject
rejection.ndf<-rest.valdes.ndf[[2]]
ave.rejection.ndf[i,]<-apply(rejection.ndf,2,mean)
# mean and 100(1-alpha)% interval of logOR
lor<-param.ndf[,7] # extract logOR
result[i,1:3] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))


param.nndf<-rest.valdes.nndf[[1]] 
ave.param.nndf[i,]<-apply(param.nndf,2,mean)
# rejection vectors: 0=accept, 1=reject
rejection.nndf<-rest.valdes.nndf[[2]]
ave.rejection.nndf[i,]<-apply(rejection.nndf,2,mean)
# mean and 100(1-alpha)% interval of logOR
lor<-param.nndf[,7] # extract logOR
result[i,4:6] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))


param.df<-rest.valdes.df[[1]] 
ave.param.df[i,]<-apply(param.df,2,mean)
# rejection vectors: 0=accept, 1=reject
rejection.df<-rest.valdes.df[[2]]
ave.rejection.df[i,]<-apply(rejection.df,2,mean)
# mean and 100(1-alpha)% interval of logOR
lor<-param.df[,7] # extract logOR
result[i,7:9] <- c(mean(lor), quantile(lor, c(alpha/2,1-alpha/2)))


} #end for i


flag1 <- ifelse((Logodd>=result[,2] & Logodd<=result[,3]),1,0)
flag2 <- ifelse((Logodd>=result[,5] & Logodd<=result[,6]),1,0)
flag3 <- ifelse((Logodd>=result[,8] & Logodd<=result[,9]),1,0)


# MSE
sumr[1,]<-c(mean((result[,1]-Logodd)^2),mean((result[,4]-Logodd)^2),mean((result[,7]-Logodd)^2))
# Coverage proportion
sumr[2,]<-c(sum(flag1)/NREP,sum(flag2)/NREP,sum(flag3)/NREP)
# Average length of the 100(1-alpha)% interval
sumr[3,]<-c(mean(result[,3]-result[,2]),mean(result[,6]-result[,5]),mean(result[,9]-result[,8]))


summary[(3*(SCEN-1)+1):(3*SCEN),]<- sumr
print(result)

if (SCEN==1) {
        write.table(result,file="Bayes/All/result1.txt",sep=" ")
        write.table(ave.param.ndf,file="Bayes/All/NDF/avepar1.txt",sep=" ")
        write.table(ave.rejection.ndf,file="Bayes/All/NDF/averejec1.txt",sep=" ")
        # save param for the NREPth dataset in each scenario (i.e. the 400th dataset in each scenario, by default)
    write.table(param.ndf,file="Bayes/All/NDF/lastpar1.txt",sep=" ")

        write.table(ave.param.nndf,file="Bayes/All/NNDF/avepar1.txt",sep=" ")
        write.table(ave.rejection.nndf,file="Bayes/All/NNDF/averejec1.txt",sep=" ")
        # save param at i=400
    write.table(param.nndf,file="Bayes/All/NNDF/lastpar1.txt",sep=" ")

        write.table(ave.param.df,file="Bayes/All/DF/avepar1.txt",sep=" ")
        write.table(ave.rejection.df,file="Bayes/All/DF/averejec1.txt",sep=" ")
        # save param at i=400
    write.table(param.df,file="Bayes/All/DF/lastpar1.txt",sep=" ")

}

if (SCEN==2) {
        write.table(result,file="Bayes/All/result2.txt",sep=" ")
        write.table(ave.param.ndf,file="Bayes/All/NDF/avepar2.txt",sep=" ")
        write.table(ave.rejection.ndf,file="Bayes/All/NDF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param.ndf,file="Bayes/All/NDF/lastpar2.txt",sep=" ")

        write.table(ave.param.nndf,file="Bayes/All/NNDF/avepar2.txt",sep=" ")
        write.table(ave.rejection.nndf,file="Bayes/All/NNDF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param.nndf,file="Bayes/All/NNDF/lastpar2.txt",sep=" ")

        write.table(ave.param.df,file="Bayes/All/DF/avepar2.txt",sep=" ")
        write.table(ave.rejection.df,file="Bayes/All/DF/averejec2.txt",sep=" ")
        # save param at i=400
    write.table(param.df,file="Bayes/All/DF/lastpar2.txt",sep=" ")
}

if (SCEN==3) {
        write.table(result,file="Bayes/All/result3.txt",sep=" ")
        write.table(ave.param.ndf,file="Bayes/All/NDF/avepar3.txt",sep=" ")
        write.table(ave.rejection.ndf,file="Bayes/All/NDF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param.ndf,file="Bayes/All/NDF/lastpar3.txt",sep=" ")

        write.table(ave.param.nndf,file="Bayes/All/NNDF/avepar3.txt",sep=" ")
        write.table(ave.rejection.nndf,file="Bayes/All/NNDF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param.nndf,file="Bayes/All/NNDF/lastpar3.txt",sep=" ")

        write.table(ave.param.df,file="Bayes/All/DF/avepar3.txt",sep=" ")
        write.table(ave.rejection.df,file="Bayes/All/DF/averejec3.txt",sep=" ")
        # save param at i=400
    write.table(param.df,file="Bayes/All/DF/lastpar3.txt",sep=" ")
}

if (SCEN==4) {
        write.table(result,file="Bayes/All/result4.txt",sep=" ")
        write.table(ave.param.ndf,file="Bayes/All/NDF/avepar4.txt",sep=" ")
        write.table(ave.rejection.ndf,file="Bayes/All/NDF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param.ndf,file="Bayes/All/NDF/lastpar4.txt",sep=" ")

        write.table(ave.param.nndf,file="Bayes/All/NNDF/avepar4.txt",sep=" ")
        write.table(ave.rejection.nndf,file="Bayes/All/NNDF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param.nndf,file="Bayes/All/NNDF/lastpar4.txt",sep=" ")

        write.table(ave.param.df,file="Bayes/All/DF/avepar4.txt",sep=" ")
        write.table(ave.rejection.df,file="Bayes/All/DF/averejec4.txt",sep=" ")
        # save param at i=400
    write.table(param.df,file="Bayes/All/DF/lastpar4.txt",sep=" ")
}


} #end for SCEN

print(summary)

write.table(summary,file="Bayes/All/summary.txt",sep=" ")

} #end misc=All


} #END OF BayesMisc