Bayesian Coverage Probability

The following code reproduces Figures of coverage probability

• Figure 5 for the OUOUCIR in the manuscript.

• Figure 2 for the OUBMCIR model in the supplemental material.

rm(list=ls())
library(bayestestR)
library(ggplot2)
library(gridExtra)

#if the url does not work, please download the RData file 
load(url("http://www.tonyjhwueng.info/ououcir/ouxxcircoverageprob.RData"))

# due to ggplot2, please run a couple time to get the figure
# OUOUCIR model
grid.arrange(p.ououcir.alpha.y,p.ououcir.alpha.x,p.ououcir.theta.x,p.ououcir.sigmasq.x, p.ououcir.alpha.tau, p.ououcir.theta.tau, p.ououcir.sigmasq.tau,p.ououcir.b0,p.ououcir.b1,p.ououcir.b2,nrow=2)

#OUBMCIR model
grid.arrange(p.oubmcir.alpha.y, p.oubmcir.sigmasq.x, p.oubmcir.alpha.tau, p.oubmcir.theta.tau, p.oubmcir.sigmasq.tau,p.oubmcir.b0,p.oubmcir.b1,p.oubmcir.b2,nrow=2)

The following summarize the models result

##################################
##########SUMMARY HERE############
##################################
setwd("http://www.tonyjhwueng.info/ououcir/simulation_32_bciV3/")
rm(list=ls())
library(bayestestR)


true.alpha.y<-0.2
true.alpha.x<-0.125
true.theta.x<-0
true.sigmasq.x<-2
true.alpha.tau<-0.25
true.theta.tau<-1.5
true.sigmasq.tau<-1
true.tau<-1

true.b0 <- 0
true.b1 <- 1
true.b2 <- 0.5


TrueParamArray<-c(
  true.alpha.y,
  true.alpha.x,
  true.theta.x,
  true.sigmasq.x,
  true.alpha.tau,
  true.theta.tau,
  true.sigmasq.tau,
  true.tau,
  true.b0,
  true.b1,
  true.b2)

names(TrueParamArray)<-c(
"alpha.y",
"alpha.x",
"theta.x",
"sigmasq.x",
"alpha.tau",
"theta.tau",
"sigmasq.tau",
"tau",
"b0",
"b1",
"b2")

modelset<-c("ououcir","oubmcir","ououbm","oubmbm")
#modelset<-c("ououcir")
percentBCI<-seq(5,100,by=5)
mainfolder<-"~/Dropbox/JournalSubmission/compstat-oubmcir/simulation_32_bciV3/"
sizearray<-c(32,64)
subfolder<-c("uniform32wider/","uniform64wider/")
#subfolder<-c("uniform32moderateB/","uniform64moderate/")
#subfolder<-c("uniform32moderateBmorenarrow/","uniform64morenarrow/")


for(modelIndex in 1:length(modelset)){
  #  modelIndex<-4
  for(subfolderIndex in 1:length(subfolder)){
    #subfolderIndex<-2
  foldertoload<-paste(mainfolder,subfolder[subfolderIndex],sep="")
  model<-modelset[modelIndex]
  load(paste(foldertoload,model,"simRep",sizearray[subfolderIndex],".RData",sep=""))
  param.inci.freq.array<-array(0,c(length(percentBCI),dim(postsample)[2]))
  head(param.inci.freq.array)
  colnames(param.inci.freq.array)<-colnames(postsample)
  rownames(param.inci.freq.array)<-paste(percentBCI,"%",sep="")
  for(repIndex in 1:reps){
#   repIndex<-1
    postsample<-get(paste(model,"Rep",repIndex,sep=""))
    #head(postsample)
#    quanbound<-apply(postsample,2,ci,ci=percentBCI,method="HDI")
#    quanbound<-as.data.frame(quanbound)
#    head(quanbound)
    for(paramname in colnames(postsample)){
#      paramname<-"alpha.y"
#      print(paramname)
      true.param<-TrueParamArray[paramname]
      for(percentIndex in 1:length(percentBCI)){
#         percentIndex<-1
         bdd<-ci(postsample[,paramname],ci=percentBCI[percentIndex]/100,method="HDI")
         lhb<-bdd$CI_low
         rhb<-bdd$CI_high
         if(true.param<rhb & true.param>lhb){
         param.inci.freq.array[percentIndex,paramname]<-param.inci.freq.array[percentIndex,paramname]+1
         }
      }
    }
  }
  assign(paste(model,"InCIfreq",sizearray[subfolderIndex],sep=""),param.inci.freq.array)
  rm(param.inci.freq.array)
}
}
  
paste(modelset,"InCIfreq",sep="")
ououcirInCIfreq<-rbind(ououcirInCIfreq32,ououcirInCIfreq64)
oubmcirInCIfreq<-rbind(oubmcirInCIfreq32,oubmcirInCIfreq64)

#ououbmInCIfreq<-rbind(ououbmInCIfreq32,ououbmInCIfreq64)
#oubmbmInCIfreq<-rbind(oubmbmInCIfreq32,oubmbmInCIfreq64)

library(ggplot2)


percentFreqbci<-function(param.inci.freq=param.inci.freq,paramname=paramname){
  percent<-seq(5,100,by=5)
  taxa<- c(rep(32,length(percent)),rep(64,length(percent)))
  percent<-rep(percent,2)
  df<-data.frame(percent=percent,taxa=taxa,param.inci.freq=param.inci.freq)
  p<- ggplot(data=df,aes(x=percent,y=param.inci.freq)) + geom_point(aes(shape=factor(taxa)),size=3.6) + geom_abline(intercept=0,slope=1,color="black",size=1.5) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + scale_x_continuous("Expected Probability", labels = as.character(percent), breaks = percent)+ scale_y_continuous("Estimated Probability", labels = as.character(percent), breaks = percent) + ggtitle(paramname)+  theme(plot.title = element_text(size = 20, face = "bold")) +
    theme(plot.title = element_text(hjust = 0.5))+ theme(legend.position = "none") #geom_bar(stat="identity")
  return(p)
}



library(gridExtra)

##############
### OUOUCIR###
##############

df <- data.frame()
p.ououcir<-ggplot(df)+ xlim(0,10)+ylim(0,10) +geom_blank()+annotate("text", x=5, y=8, label= "Model",size=8)+annotate("text", x=5, y=6, label= "OUOUCIR",size=8)+theme_void()

param.inci.freq<-ououcirInCIfreq[,"alpha.y"]
paramname<-expression(paste( alpha[y],sep= "     "))
p.ououcir.alpha.y<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-ououcirInCIfreq[,"alpha.x"]
paramname<-expression(paste( alpha[x],sep= "     "))
p.ououcir.alpha.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"theta.x"]
paramname<-expression(paste( theta[x],sep= "     "))
p.ououcir.theta.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"alpha.tau"]
paramname<-expression(paste( alpha[tau],sep= "     "))
p.ououcir.alpha.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"theta.tau"]
paramname<-expression(paste( theta[tau],sep= "     "))
p.ououcir.theta.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"sigmasq.tau"]
paramname<-expression(paste( sigma[tau]^2,sep= "     "))
p.ououcir.sigmasq.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"sigmasq.x"]
paramname<-expression(paste( sigma[x]^2,sep= "     "))
p.ououcir.sigmasq.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"b0"]
paramname<-expression(paste(beta[0],sep= "     "))
p.ououcir.b0<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"b1"]
paramname<-expression(paste(beta[1],sep= "     "))
p.ououcir.b1<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououcirInCIfreq[,"b2"]
paramname<-expression(paste(beta[2],sep= "     "))
p.ououcir.b2<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)



grid.arrange(p.ououcir.alpha.y,p.ououcir.alpha.x,p.ououcir.theta.x,p.ououcir.sigmasq.x, p.ououcir.alpha.tau, p.ououcir.theta.tau, p.ououcir.sigmasq.tau,p.ououcir.b0,p.ououcir.b1,p.ououcir.b2,nrow=2)




# save(
#   ououcirInCIfreq,
#   oubmcirInCIfreq,
#   p.ououcir,
#   p.ououcir.alpha.y,
#   p.ououcir.alpha.x,
#   p.ououcir.theta.x,
#   p.ououcir.alpha.tau,
#   p.ououcir.theta.tau,
#   p.ououcir.sigmasq.tau,
#   p.ououcir.sigmasq.x,
#   p.ououcir.b0,
#   p.ououcir.b1,
#   p.ououcir.b2,
#   p.oubmcir,
#   p.oubmcir.alpha.y,
#   p.oubmcir.alpha.tau,
#   p.oubmcir.theta.tau,
#   p.oubmcir.sigmasq.tau,
#   p.oubmcir.sigmasq.x,
#   p.oubmcir.b0,
#   p.oubmcir.b1,
#   p.oubmcir.b2,file="ouxxcircoverageprob.RData")




##############
### OUBMCIR###
##############

df <- data.frame()
p.oubmcir<-ggplot(df)+ xlim(0,10)+ylim(0,10) +geom_blank()+annotate("text", x=5, y=8, label= "Model",size=8)+annotate("text", x=5, y=6, label= "OUBMCIR",size=8)+theme_void()

param.inci.freq<-oubmcirInCIfreq[,"alpha.y"]
paramname<-expression(paste( alpha[y],sep= "     "))
p.oubmcir.alpha.y<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmcirInCIfreq[,"alpha.tau"]
paramname<-expression(paste( alpha[tau],sep= "     "))
p.oubmcir.alpha.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmcirInCIfreq[,"theta.tau"]
paramname<-expression(paste( theta[tau],sep= "     "))
p.oubmcir.theta.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmcirInCIfreq[,"sigmasq.tau"]
paramname<-expression(paste( sigma[tau]^2,sep= "     "))
p.oubmcir.sigmasq.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmcirInCIfreq[,"sigmasq.x"]
paramname<-expression(paste( sigma[x]^2,sep= "     "))
p.oubmcir.sigmasq.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmcirInCIfreq[,"b0"]
paramname<-expression(paste(beta[0],sep= "     "))
p.oubmcir.b0<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmcirInCIfreq[,"b1"]
paramname<-expression(paste(beta[1],sep= "     "))
p.oubmcir.b1<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmcirInCIfreq[,"b2"]
paramname<-expression(paste(beta[2],sep= "     "))
p.oubmcir.b2<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)



grid.arrange(p.oubmcir.alpha.y, p.oubmcir.sigmasq.x, p.oubmcir.alpha.tau, p.oubmcir.theta.tau, p.oubmcir.sigmasq.tau,p.oubmcir.b0,p.oubmcir.b1,p.oubmcir.b2,nrow=2)



##############
### OUOUBM ###
##############

df <- data.frame()
p.ououbm<-ggplot(df)+ xlim(0,10)+ylim(0,10) +geom_blank()+annotate("text", x=5, y=8, label= "Model",size=12)+annotate("text", x=5, y=6, label= "OUOUBM",size=10)+theme_void()

param.inci.freq<-ououbmInCIfreq[,"alpha.y"]
paramname<-expression(paste( alpha[y],sep= "     "))
p.ououbm.alpha.y<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-ououbmInCIfreq[,"alpha.x"]
paramname<-expression(paste( alpha[x],sep= "     "))
p.ououbm.alpha.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-ououbmInCIfreq[,"theta.x"]
paramname<-expression(paste( theta[x],sep= "     "))
p.ououbm.theta.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmbmInCIfreq[,"sigmasq.x"]
paramname<-expression(paste( sigma[x]^2,sep= "     "))
p.ououbm.sigmasq.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmbmInCIfreq[,"tau"]
paramname<-expression(paste( tau,sep= "     "))
p.ououbm.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououbmInCIfreq[,"b0"]
paramname<-expression(paste(beta[0],sep= "     "))
p.ououbm.b0<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououbmInCIfreq[,"b1"]
paramname<-expression(paste(beta[1],sep= "     "))
p.ououbm.b1<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-ououbmInCIfreq[,"b2"]
paramname<-expression(paste(beta[2],sep= "     "))
p.ououbm.b2<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


##############
### OUBMBM ###
##############

df <- data.frame()
p.oubmbm<-ggplot(df)+ xlim(0,10)+ylim(0,10) +geom_blank()+annotate("text", x=5, y=8, label= "Model",size=12)+annotate("text", x=5, y=6, label= "OUBMBM",size=10)+theme_void()

param.inci.freq<-oubmbmInCIfreq[,"alpha.y"]
paramname<-expression(paste( alpha[y],sep= "     "))
p.oubmbm.alpha.y<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmbmInCIfreq[,"sigmasq.x"]
paramname<-expression(paste( sigma[x]^2,sep= "     "))
p.oubmbm.sigmasq.x<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)

param.inci.freq<-oubmbmInCIfreq[,"tau"]
paramname<-expression(paste( tau,sep= "     "))
p.oubmbm.tau<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmbmInCIfreq[,"b0"]
paramname<-expression(paste(beta[0],sep= "     "))
p.oubmbm.b0<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmbmInCIfreq[,"b1"]
paramname<-expression(paste(beta[1],sep= "     "))
p.oubmbm.b1<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


param.inci.freq<-oubmbmInCIfreq[,"b2"]
paramname<-expression(paste(beta[2],sep= "     "))
p.oubmbm.b2<-percentFreqbci(param.inci.freq=param.inci.freq,paramname=paramname)


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


grid.arrange(p.ououcir.alpha.y,p.ououcir.alpha.x,p.ououcir.theta.x,p.ououcir.sigmasq.x, p.ououcir.alpha.tau, p.ououcir.theta.tau, p.ououcir.sigmasq.tau,p.ououcir.b0,p.ououcir.b1,p.ououcir.b2,nrow=2)

grid.arrange(p.oubmcir.alpha.y, p.oubmcir.sigmasq.x, p.oubmcir.alpha.tau, p.oubmcir.theta.tau, p.oubmcir.sigmasq.tau,p.oubmcir.b0,p.oubmcir.b1,p.oubmcir.b2,nrow=2)

grid.arrange(p.ououbm.alpha.y, p.ououbm.alpha.x, p.ououbm.theta.x, p.ououbm.sigmasq.x, p.ououbm.tau,p.ououbm.b0,p.ououbm.b1,p.ououbm.b2,nrow=2)

grid.arrange(p.oubmbm.alpha.y,p.oubmbm.sigmasq.x,p.oubmbm.tau,p.oubmbm.b0,p.oubmbm.b1,p.oubmbm.b2,nrow=2)