## ----warning=FALSE, message=FALSE------------------------------------------
library("ggplot2")
library("dplyr")
library("tidyr")
library("gridExtra")
#library("LSD")
library("IHW")

## --------------------------------------------------------------------------
# try function with non-monotonic relation
sim <- function(m =20000){
  x <- 1-runif(m)
  pi0 <- 0.75+0.3*(x-0.5)^2#0.8 #ifelse(x<=.5, 0.5 + 0.6*x, 0.8)
  eff_size <- 2-2.2*abs(x-0.5)^1.3  #was 3.2 #ifelse(x <=.5, 1.5, 1.5-(x-1/2)*4)
  H <- rbinom(m, 1, 1-pi0)
  Z <- rnorm(m)
  Z[H==1] <- rnorm(sum(H), eff_size[H==1])
  pvalue <- 1-pnorm(Z)
  t <- pvalue
  t_bh <- get_bh_threshold(pvalue,0.1)

  alt_ft_bh <- dnorm(-eff_size+ qnorm(1-t_bh))/dnorm(qnorm(1-t_bh))
  tdr_bh <- alt_ft_bh*(1-pi0)/((1-pi0)*alt_ft_bh + pi0)
  t_ddhw <- thresholds(ihw(pvalue, x, alpha= .1, nbins=20,
                           nsplits_internal=10, nfolds = 1), levels_only=FALSE)

  alt_ft_ddhw <- dnorm(-eff_size+ qnorm(1-t_ddhw))/dnorm(qnorm(1-t_ddhw))
  tdr_ddhw <- alt_ft_ddhw*(1-pi0)/((1-pi0)*alt_ft_ddhw + pi0)

  alt_ft <- dnorm(qnorm(1-t)-eff_size)/dnorm(qnorm(1-t))
  tdr <- alt_ft*(1-pi0)/((1-pi0)*alt_ft + pi0)
  return(data.frame(x=x,pi0=pi0, t=t, eff_size=eff_size,H=H,
                    Z=Z, pvalue=pvalue, alt_ft=alt_ft, tdr=tdr,
                    tdr_bh=tdr_bh, alt_ft_bh=alt_ft_bh, t_bh=t_bh,
                    t_ddhw=t_ddhw, tdr_ddhw=tdr_ddhw, alt_ft_ddhw=alt_ft_ddhw))
}

## --------------------------------------------------------------------------
set.seed(1)
sim_df <- sim(m=80000)

## --------------------------------------------------------------------------
sim_df_t <- select(sim_df, x, t_bh, t_ddhw ) %>%
            gather(method, t, -x) %>%
            mutate(method=ifelse(method=="t_bh", "BH","IHW"))

ggplot(sim_df_t, aes(x=x, y=t,color=method)) +
                   geom_step() +
                   xlab("covariate") +
                   ylab("p-value threshold")

## --------------------------------------------------------------------------
sim_df_tdr <- select(sim_df, x, tdr_bh, tdr_ddhw ) %>%
              gather(method, tdr, -x) %>%
              mutate(method=ifelse(method=="tdr_bh", "BH","IHW"))


ggplot(sim_df_tdr, aes(x=x, y=tdr,color=method)) +
                   geom_step() +
                   xlab("covariate") +
                   ylab("tdr threshold")

## --------------------------------------------------------------------------
#sim_df$group <- groups_by_filter(sim_df$x, 10)
#plot(sim_df$x, sim_df$tdr_bh)

#with(filter(sim_df,-log10(pvalue) > 0.001), heatscatter(-log10(pvalue),x, add.contour=TRUE))
#with(filter(sim_df,tdr>0.2), heatscatter(tdr,x, add.contour=TRUE)) 

#ggplot(filter(sim_df), aes(x=-log10(pvalue), y=x)) + geom_density2d(aes(color = ..level..)) 
#ggplot(filter(sim_df, tdr>0.4), aes(x=tdr, y=x)) + geom_density2d(aes(color = ..level..))