##-----------------------------------------------------------------------------
## ST790 -- R codes for Lecture 10a
## Author: Ryan Martin (rgmarti3@ncsu.edu)                                      
##-----------------------------------------------------------------------------

# Marginal (partial-prior) IM for the odds ratio in a 2x2 table

or.pmf <- function(th, x, y, n) {
  
  n1 <- n[1]
  n2 <- n[2]
  lo <- max(x - n1, 0)
  hi <- min(n2, x)
  s <- lo:hi
  num <- choose(n2, s) * choose(n1, x - s) * exp(s * th)
  den <- sum(num)
  return(num[which(lo:hi == y)] / den)
  
}

or.log.rl <- function(th, x, y, n) {
  
  th1 <- y / n[1]
  th2 <- (x - y) / n[2]
  hat <- 0 #log(th1) - log(1 - th1) - log(th2) + log(1 - th2)
  f <- function(z) -log(or.pmf(z, x, y, n))
  opt <- nlm(f, hat)$minimum
  return(-f(th) + opt)
  
}

or.pl.vac <- function(th, x, y, n) {
  
  n1 <- n[1]
  n2 <- n[2]
  lo <- max(x - n1, 0)
  hi <- min(n2, x)
  s <- lo:hi
  z <- or.log.rl(th, x, y, n)
  Z <- sapply(s, function(m) or.log.rl(th, x, m, n))
  p <- sapply(s, function(m) or.pmf(th, x, m, n))
  return(sum(p[Z <= z]))
  
}


or.pl.part <- function(th, x, y, n, qq) {
  
  n1 <- n[1]
  n2 <- n[2]
  lo <- max(x - n1, 0)
  hi <- min(n2, x)
  f <- function(z, u) log(or.pmf(u, x, z, n)) + log(qq(u))
  s <- lo:hi
  opt <- 0 * s
  for(i in seq_along(s)) {
    
    g <- function(u) -f(s[i], u)
    opt[i] <- -nlm(g, 0)$minimum
    
  }
  eta.y <- f(y, th) - opt[which(s == y)]
  h <- function(u) {
    
    eta.Y <- sapply(s, function(b) f(b, u)) - opt
    pmf <- sapply(s, function(m) or.pmf(u, x, m, n))
    return(sum(pmf[eta.Y <= eta.y]))
    
  }
  tt <- seq(-5, 5, length=201)
  pl <- sapply(tt, h)
  pr <- qq(tt)
  k <- function(a) sapply(a, function(z) max(pr[pl >= z]))
  oo <- min(pl) + integrate(k, min(pl), max(pl))$value
  return(oo)
  
}


y <- c(1, 2)
n <- c(43, 39)
th <- seq(-4, 7, length=100)
qq <- function(u) pmin(1, 1 / abs(u))
pl.vac <- sapply(th, function(v) or.pl.vac(v, sum(y), y[2], n))
pl.part <- sapply(th, function(v) or.pl.part(v, sum(y), y[2], n, qq))
plot(th, pl.vac, type="l", col="white", xlab=expression(theta), ylab="Plausibility")
curve(qq, col=2, lty=3, add=TRUE)
lines(th, pl.part, col=2)
lines(th, pl.vac)