tests/testthat/test_im04ci.R

In matteobonvini/sensitivitypuc: Sensitivity Analysis via the Proportion of Unmeasured Confounding (PUC)

context("Confidence Interval from Imbens & Manski (2004)")

test_that("Coverage of Imbens & Manski (2004) is correct", {
  
  require(pbapply)
  # Data generating process:
  # X ~ Unif(0, 1); A ~ Bern(0.5); Y = A*X + (1-A)
  # so that pi(x) = 0.5; mu1(x) = x; mu0(x) = 1; g(x) = 1 - 0.5*x
  # eps-quantile of g(x) = 0.5*(1+eps)
  # E(mu1(x) - mu0(x)) = -0.5
  # upper bound is psiu(eps) = -0.25(eps^2 - 4*eps + 2)
  # hence eps0 is 0.5*(4 - sqrt(8)) = 0.5858
  # lower bound is psil(eps) = 0.25(eps^2 - 2*eps - 2)
  
  eps <- seq(0, 0.1, 0.001)
  
  # Start simulation to check coverage
  nsim <- 1000
  n <- 1000
  alpha <- 0.2
  
  psiu <- -0.25*(eps^2 - 4*eps + 2)
  psil <- 0.25*(eps^2 - 2*eps - 2)
  
  sim_fn <- function() {
    # Simulation function to estimate eps0
    # Generate data according to the model above
    x <- runif(n, 0, 1)
    a <- rbinom(n, 1, 0.5)
    y <- a*x + (1-a)
    
    # Regression functions are known exacty
    mu1x <- x
    mu0x <- rep(1, n)
    pi1x <- pi0x <- rep(0.5, n)
    nuis_fns <- matrix(NA, ncol = 7, nrow = 2 * n,
                       dimnames = list(NULL, c("mu1", "mu0", "pi1", "pi0",
                                               "unit_num", "fold", "is_test")))
    nuis_fns[, "mu1"] <- rep(mu1x, 2)
    nuis_fns[, "mu0"] <- rep(mu0x, 2)
    nuis_fns[, "pi1"] <- rep(pi1x, 2)
    nuis_fns[, "pi0"] <- rep(pi0x, 2)
    
    nuis_fns[, "unit_num"] <- rep(1:n, 2)
    nuis_fns[, "fold"] <- 1
    nuis_fns[, "is_test"] <- c(rep(1, n), rep(0, n))
    
    nuis_fns_list <- list(test = nuis_fns[1:n, ], 
                          train = nuis_fns[(n+1):(2*n), ],
                          order_obs = nuis_fns[1:n, "unit_num"])
    gx <- 1 - 0.5*x
    
    res <- get_bound(y = y, a = a, x = as.data.frame(x), ymin = 0, ymax = 1, 
                     outfam = NULL, treatfam = NULL, model = "x", eps = eps, 
                     delta = 1, nsplits = 1, do_mult_boot = FALSE,
                     do_eps_zero = FALSE, nuis_fns = nuis_fns_list, alpha = alpha,
                     do_rearrange = FALSE)
    # Machine precision 
    
    calpha <- round(res$im04_calpha, 15)
    zalpha <- round(qnorm(1-alpha), 15)
    zalpha2 <- round(qnorm(1-alpha/2), 15)
    expect_true(all(zalpha <= calpha & calpha <= zalpha2))
    
    ci_lo_im04 <- round(res$bounds[, "ci_im04_lo", 1], 15)
    ci_hi_im04 <- round(res$bounds[, "ci_im04_hi", 1], 15)
    ci_lo_pt <- round(res$bounds[, "ci_lb_lo_pt", 1], 15)
    ci_hi_pt <- round(res$bounds[, "ci_ub_hi_pt", 1], 15)
    
    # IM 04 Ci should be contained in the pointwise bands covering the 
    # identification region
    expect_true(all(ci_lo_pt <= ci_lo_im04 & ci_hi_im04 <= ci_hi_pt))
    
    # prepare output
    out <- matrix(c(ci_lo_pt, ci_hi_pt, ci_lo_im04, ci_hi_im04, calpha),
                  ncol = 5, nrow = length(eps),
                  dimnames = list(NULL, c("ci_lo_pt", "ci_lo_pt", "ci_lo_im04", 
                                          "ci_hi_im04", "calpha")))
    return(out)
  }
  
  cvg <- Vectorize(function(sims, psi, x) {
    coverage(sims[x, "ci_lo_im04", ], sims[x, "ci_hi_im04", ], 
             psi[x], psi[x])
  }, vectorize.args = "x")
  
  # Simulation begins
  sims <- pbreplicate(nsim, sim_fn())
  
  # Check pointwise coverage is okay at, e.g., eps = eps[5]
  # case when the true curve psi(eps) is equal to lower bound
  cvgs_l <- cvg(sims = sims, psi = psil, x = 1:length(eps))
  # case when the true curve psi(eps) is equal the upper bound
  cvgs_u <- cvg(sims = sims, psi = psiu, x = 1:length(eps))
  # Both coverage and bias are expressed in %
  expect_true(all(abs(cvgs_l - 100*(1-alpha)) <= 5))
  expect_true(all(abs(cvgs_u - 100*(1-alpha)) <= 5))
})