R/mWald_p.R

In BriceonWiley/IntegratedLikelihood.R:

mWald_p <- function(n00 = 49, n01 = 1, n10 = 1, n11 = 2, X = 14, Y = 433,
                    alpha = 0.05) {

  sample <- c(n00, n01, n10, n11, X, Y)
  mles <- calculate_mles_binom(sample)
  pihat <- mles$pihat
  lam1hat <- mles$lam1hat
  lam2hat <- mles$lam2hat
  phat <- mles$phat
  s_hat <- log(phat / (1 - phat))

  frac1 <- (pihat * lam1hat * (1 - lam1hat)) / sum(sample[1:4])
  frac2 <- (1 - pihat) * lam2hat * (1 - lam2hat) / sum(sample[1:4])
  frac3 <- (lam1hat - lam2hat)^2 * pihat * (1 - pihat) / sum(sample)
  sigmahat2 <- frac1 + frac2 + frac3
  tauhat <- sigmahat2 / (phat * (1 - phat))^2

  zcrit <- qnorm(1 - alpha / 2)
  L <- exp(s_hat - zcrit * sqrt(tauhat)) / (exp(s_hat - zcrit * sqrt(tauhat)) + 1)
  U <- exp(s_hat + zcrit * sqrt(tauhat)) / (exp(s_hat + zcrit * sqrt(tauhat)) + 1)
  tibble("mW_Lower" = L, "mW_Upper" = U)
}