R/coverage_gam.R

In BriceonWiley/IntegratedLikelihood.R:

coverage_gam <- function(x, ignore_NaN = FALSE) {
  N1 <- x$counts["N1"]
  N2 <- x$counts["N2"]
  p1 <- x$proportions["p1"]
  p2 <- x$proportions["p2"]
  NaNs <- x$counts["NaN"]
  gam <- log((p1 * (1 - p2)) / (p2 * (1 - p1)))
  int_names <- names(x$intervals)
  if (ignore_NaN) {
    # usable <- filter(x$intervals, !is.nan(W_Lower)) %>%
    #   filter(!is.nan(ILR_Lower))                        # do this if ILR is NA
    usable <- filter(x$intervals, !is.nan(W_Lower))
    lower <- select(usable, str_subset(int_names, "Lower"))
    upper <- select(usable, str_subset(int_names, "Upper"))
  } else {
    lower <- select(x$intervals, str_subset(int_names, "Lower"))
    upper <- select(x$intervals, str_subset(int_names, "Upper"))
  }
  coverages <- map2(lower, upper, ~ map2(.x, .y, ~ between(gam, .x, .y))) %>%
    transpose() %>%
    bind_rows() %>%
    summarise_all(mean)
  names(coverages) <- str_remove(names(coverages), "[_]?Lower[_]?")
  if (N1 == N2) {
    coverages %>%
      add_column(
        "N" = N1, "gamma" = gam, "p1" = p1, "p2" = p2, "NaNs" = NaNs,
        .before = TRUE
      )
  } else {
    coverages %>%
      add_column(
        "N1" = N1, "N2" = N2,
        "gamma" = gam, "p1" = p1, "p2" = p2, "NaNs" = NaNs,
        .before = TRUE
      )
  }
}