Qbeta_Issues.R

library(tidyverse)
library(ggplot2); theme_set(theme_bw())
library(viridis)
library(bbmle)
library(DPQ)
library(deSolve)
source("testing_funs.R")

### Testing Qbeta2 vs qbeta
dd <- (expand.grid(test_prop=c(0.001,0.0012,0.0015,0.002,0.005,0.01,0.05,0.1,0.25),
                   phi=seq(from=0.01,to=0.99,by=0.01),
                   inc=seq(from=0.01,to=0.99,by=0.01))
       %>% as_tibble()
       %>% mutate(pos_prop_qlog=prop_pos_test_new(inc,test_prop,phi,method="est",qfun=qbeta))
       %>% mutate(pos_prop_Qlog=prop_pos_test_new(inc,test_prop,phi,method="est"))
       %>% mutate(diff_Qq_dif=abs_log_sub(pos_prop_qlog,pos_prop_Qlog))
       )

#which(dd$diff_Qq_dif!=-Inf)

fig_logdiff_Qbeta_qbeta <- (
  ggplot(dd,aes(phi,inc,fill=diff_Qq_dif,group=test_prop))
  + geom_raster()
  + facet_wrap(~test_prop,scale="free",labeller = label_both)
  #+ scale_y_log10()
  + scale_fill_viridis_c(expand=c(0,0))
  + ggtitle("log-space diff between qbeta+est and Qbeta2+est")
)
ggsave("log-diff_Qbeta-qbeta.png",plot=fig_logdiff_Qbeta_qbeta, path = "./pix", width=3200,height=1800,units="px")

print(ggplot(dd,aes(phi,inc,fill=pos_prop_qlog,group=test_prop))
      + geom_raster()
      + facet_wrap(~test_prop,scale="free",labeller = label_both)
      #+ scale_y_log10()
      + scale_fill_viridis_c(expand=c(0,0))
      + ggtitle("expected positivity generated by qbeta+log")
)

print(ggplot(dd,aes(phi,inc,fill=pos_prop_Qlog,group=test_prop))
      + geom_raster()
      + facet_wrap(~test_prop,scale="free",labeller = label_both)
      #+ scale_y_log10()
      + scale_fill_viridis_c(expand=c(0,0))
      + ggtitle("expected positivity generated by Qbeta2+log")
)

##Seems in practice qbeta does not create more trouble than Qbeta2