R/intLik_gam.R

In BriceonWiley/IntegratedLikelihood.R:

intLik_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()


  # data preparation
  crit <- qchisq(1 - alpha, 1)
  lb <- Wald_gam(sample1, sample2)[[1]]
  ub <- Wald_gam(sample1, sample2)[[2]]


  # check if CI is possible
  if (any(is.na(c(lb, ub)))) return(tibble("ILR_Lower" = NaN, "ILR_Upper" = NaN))


  # find MLE
  gamhat <- optimize(
    likelihood_gam, c(lb - 2, ub + 2), data1 = sample1, data2 = sample2,
    maximum = TRUE
    )$maximum
  num <- likelihood_gam(gamhat, sample1, sample2)
  lb <- gamhat - 2 * crit; ub <- gamhat + 2 * crit


  # build CI
  lower <- tryCatch(
    uniroot(
      likelihood_ratio_gam,
      lower = lb,
      upper = gamhat,
      num = num,
      data1 = sample1,
      data2 = sample2,
      crit = crit
      )$root,
    error = function(e) {
      message(
        paste0(
          "There was a problem with the lower bound for the sample:\n",
          "Sample 1: n00 = ", sample1[1], ", n01 = ", sample1[2],
          ", n10 = ", sample1[3], ", n11 = ", sample1[4], ", X = ", sample1[5],
          " and Y = ", sample1[6], ".\n", "Sample 2: n00 = ", sample2[1],
          ", n01 = ", sample2[2], ", n10 = ", sample2[3], ", n11 = ",
          sample2[4], ", X = ", sample2[5], " and Y = ", sample2[6], "."
        )
      )
      NaN
    }
    )
  upper <- tryCatch(
    uniroot(
      likelihood_ratio_gam,
      lower = gamhat,
      upper = ub,
      num = num,
      data1 = sample1,
      data2 = sample2,
      crit = crit
      )$root,
    error = function(e) {
      message(
        paste0(
          "There was a problem with the upper bound for the sample:\n",
          "Sample 1: n00 = ", sample1[1], ", n01 = ", sample1[2],
          ", n10 = ", sample1[3], ", n11 = ", sample1[4], ", X = ", sample1[5],
          " and Y = ", sample1[6], ".\n", "Sample 2: n00 = ", sample2[1],
          ", n01 = ", sample2[2], ", n10 = ", sample2[3], ", n11 = ",
          sample2[4], ", X = ", sample2[5], " and Y = ", sample2[6], "."
        )
      )
      NaN
    }
    )


  # return CI
  tibble("ILR_Lower" = lower, "ILR_Upper" = upper)
}