Bounds given different values of \(\tau\) and \(\rho\)

markov  <- function(t, r) matrix(c(t-r+1, 1-t-r, 1-t+r, 1+t+r)/2, 2, 2)

bounds_T <- function(T) {
  L       <- function(T) max(0, (T[2,2] - T[1,2])/T[2,2])
  H       <- function(T) min(1, T[1,1]/T[2,2])
  
  tau <- T[2,2] - T[1,2]
  rho <- T[2,2] - T[1,1]
  
  T2  <- markov(tau^.5, rho/(1+tau^.5))
  
  best_het <- function(T) matrix(c(
    min(1, T[1,2]*2), 
    2*T[1,2] - min(1, T[1,2]*2), 
    2*T[2,2] - min(1, T[2,2]*2),
    min(1, T[2,2]*2)), 
    2, 2)
  
  best_het_bounds <- function(T){
    M <- best_het(T)
    M1 <- matrix(c(1 - M[1,], M[1,]), 2, 2)
    M2 <- matrix(c(1 - M[2,], M[2,]), 2, 2)
    c(L=  (M[2,2]/(M[1,2] + M[2,2]))*L(M2)+
        (M[1,2]/(M[1,2] + M[2,2]))*L(M1),
      H = (M[2,2]/(M[1,2] + M[2,2]))*H(M2) + 
        (M[1,2]/(M[1,2] + M[2,2]))*H(M1))
  } 
  
  
  rbind(simple = c(L(T), H(T)),
        unobserved = c(L(T), L(T)),
        monotone = c(L(T), L(T)),
        two_step_homog = c(L(T2)^2, H(T2^2)),
        inifinte_homog = c(tau^{.5*(1 + rho/(1-tau))}, min(1,       tau^{rho/(1-tau)})),
        best_lower = c(T[1,1], T[1,1]),
        best_cov_un =  c(L = H(T), H = H(T)),
        best_cov_ob =  c(1,1)
  )
}

bound_names =  c("Simple bounds",
                 "Unobserved mediators",
                 "Monotonic",
                 "Two step homogeneous",
                 "Infinite step homogeneous",
                 "Best (positive) mediator",
                 "Best (unobserved) covariate",
                 "Best (observed) covariate")

# rho_values <-  round(c(-2/3, -1/3, 0, 1/3, 2/3),2)
rho_values <-  round(c(-.5, 0, .5),2)


bounds_list <- lapply(rho_values,
                      function(rho){
                        do.call(rbind, lapply(c(.1, .25), function(tau){
                          T <- markov(tau, rho)
                          bounds <- as.data.frame(bounds_T(T))
                          bounds$rho   <- as.factor(rho)
                          bounds$tau   <- as.factor(tau)
                          bounds$label <- rownames(bounds)
                          bounds
                        }))
                      })

df <- do.call(rbind, bounds_list)

df$names <- rep(bound_names, nrow(df)/length(unique(df$label)))

# fig() gives a few options for plotting:
# `bw = TRUE` plots grey colors
# `grid = TRUE` plots grid with values of tau (0.1 and 0.5) as columns and rho values as rows
# note : `grid = FALSE` requires subsetting data to given value of tau
# shape = TRUE shows shape aesthetics
fig <- function(dat = df){
  
  dat$names <- factor(dat$names, levels = rev(bound_names))
  
dat$rho <- paste("rho ==", dat$rho); 
dat$tau <- paste("tau ==", dat$tau)
dat$rho <- factor(dat$rho, levels = paste0("rho == ", rho_values))


  ggplot(data=dat, aes(x = names, ymin = L, ymax = H, color=names)) + 
    geom_point(aes(y = L), position = position_dodge(.8)) +
    geom_point(aes(y = H), position = position_dodge(.8)) +
    geom_linerange(position = position_dodge(.8)) +
    scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1)) +
    coord_flip(ylim = c(0,1)) +
    theme_bw() +
    guides(col = guide_legend(ncol = 2)) +
    labs(y = "Probability of causation", x = "", shape = "\u03C1", colour = "\u03C1") +
    theme(legend.position = "bottom") + 
    facet_grid(rho ~ tau, labeller = label_parsed) +  
    theme(legend.position = "none") + 
    scale_color_grey(start = .6, end = 0)
}


fig()