R/Wald_gam.R

In BriceonWiley/IntegratedLikelihood.R:

Wald_gam <- function(
  sample1 = c(n00 = 11, n01 = 1, n10 = 1, n11 = 2, X = 34, Y = 101),
  sample2 = c(n00 = 10, n01 = 2, n10 = 1, n11 = 2, X = 25, Y = 110),
  alpha = 0.05
) {
  # browser()
  mles1 <- calculate_mles_binom(sample1)
  mles2 <- calculate_mles_binom(sample2)

  var_gam <-
    list(mles1, mles2) %>%
    map_dbl(~ {
      f1 <-
        (.x$lam1hat - .x$lam2hat)^2 * .x$pihat * (1 - .x$pihat) / sum(.x$sample)
      f2 <-
        .x$pihat^2 * .x$lam1hat * (1 - .x$lam1hat) / sum(.x$sample[c(2, 4)])
      f3 <-
        (1 - .x$pihat)^2 * (1 - .x$lam2hat) * .x$lam2hat / sum(.x$sample[c(1, 3)])
      denom <- (.x$phat * (1 - .x$phat))^(-2)
      (f1 + f2 + f3) * denom
    }) %>%
    sum()

  gamhat <- log((mles1$phat * (1 - mles2$phat)) / (mles2$phat * (1 - mles1$phat)))
  me <- qnorm(1 - alpha / 2) * sqrt(var_gam)

  L <- gamhat - me; U <- gamhat + me
  tibble("W_Lower" = L, "W_Upper" = U)
}