test-all.R
In eco: Ecological Inference in 2x2 Tables

rm(list=ls())
library(eco)
library(testthat)
context("tests eco")

accuracy1 <- ifelse(capabilities("long.double"), 0.002, 0.005)
accuracy2 <- ifelse(capabilities("long.double"), 0.2, 0.5)

# set random seed
set.seed(12345)

# for the tests that may take a long time to finish, skip them
donotrun = 1

test_that("tests eco on registration data", {
  ## load the data
  data(reg)

  # fit the parametric model with the default prior specification
  res <- eco(Y ~ X, data = reg, verbose = TRUE)

  # summarize the results
  x <- summary(res)
  expect_that(length(x), is_equivalent_to(8))
  expect_true("W2.table" %in% names(x))
  expect_equal(x$param.table[2,1], 2.976173, tolerance = accuracy1)
  expect_equal(x$param.table[3,4], 8.238363, tolerance = accuracy1)

  # obtain out-of-sample prediction
  out <- predict(res, verbose = TRUE)
  # summarize the results
  x <- summary(out)
  expect_that(length(x), is_equivalent_to(2))
  expect_true("n.draws" %in% names(x))
  expect_equal(x$W.table[1,1], 0.4896912, tolerance = accuracy1)
  expect_equal(x$W.table[2,3], 0.3071461, tolerance = accuracy1)
})  
  
if (!donotrun) test_that("tests eco on Robinson census", {
  # load the Robinson's census data
  data(census)

  # fit the parametric model with contextual effects and N using the default prior specification
  res1 <- eco(Y ~ X, N = N, context = TRUE, data = census, verbose = TRUE)

  # summarize the results
  x <- summary(res1)
  expect_that(length(x), is_equivalent_to(8))
  expect_true("W2.table" %in% names(x))
  expect_equal(x$param.table[2,3], 2.068228, tolerance = accuracy1)
  expect_equal(x$agg.wtable[1,3], 0.6631372, tolerance = accuracy1)

  # obtain out-of-sample prediction
  out1 <- predict(res1, verbose = TRUE)
  # summarize the results
  x <- summary(out1)
  expect_that(length(x), is_equivalent_to(2))
  expect_true("n.draws" %in% names(x))
  expect_equal(x$W.table[1,3], 0.4499757, tolerance = accuracy1)
  expect_equal(x$W.table[3,1], 0.3400277, tolerance = accuracy1)
})

# ecoBD
test_that("tests ecoBD on registration data", {
  # load the registration data
  data(reg)
  # calculate the bounds
  x <- ecoBD(Y ~ X, N = N, data = reg)
  expect_that(length(x), is_equivalent_to(12))
  expect_true("aggWmin" %in% names(x))
  expect_equal(x$aggWmin[1,1], 0.216785, tolerance = accuracy1)
  expect_that(x$aggNmax[2,2], is_equivalent_to(2046800))
})

# ecoML
if (!donotrun) test_that("tests ecoML on census data", {
  # load the census data
  data(census)

  # fit the parametric model with the default model specifications
  res <- ecoML(Y ~ X, data = census[1:100,], N=census[1:100,3], epsilon=10^(-6), verbose = TRUE)
  # summarize the results
  x <- summary(res)
  expect_that(length(x), is_equivalent_to(13))
  expect_true("iters.sem" %in% names(x))
  expect_equal(x$loglik, -70.80674, tolerance = accuracy2)
  expect_equal(x$param.table[2,3], 0.07192878, tolerance = accuracy1)

  #####################################################################################
  # NOTE: this example does not work! There is no predict.ecoML defined in the package. 
  #
  # obtain out-of-sample prediction
  # out <- predict(res, verbose = TRUE)
  # summarize the results
  # summary(out)
  #####################################################################################
  
  # fit the parametric model with some individual 
  # level data using the default prior specification
  surv <- 1:600
  res1 <- ecoML(Y ~ X, context = TRUE, data = census[-surv,], supplement = census[surv,c(4:5,1)], maxit=100, verbose = TRUE)
  # summarize the results
  x <- summary(res1)
  expect_that(length(x), is_equivalent_to(13))
  expect_true("iters.sem" %in% names(x))
  expect_equal(x$loglik, -3481.877, tolerance = accuracy2)
  expect_equal(x$param.table[2,3], 0.006055498, tolerance = accuracy1)
  expect_true(is.na(x$param.table[2,2]))
  expect_true(is.na(x$param.table[2,5]))
  expect_false(is.na(x$param.table[2,6]))
})

# set random seed
set.seed(12345)

# ecoNP
test_that("tests ecoNP on census data", {
  # load the registration data
  data(reg)

  # NOTE: We set the number of MCMC draws to be a very small number in
  # the following examples; i.e., convergence has not been properly
  # assessed. See Imai, Lu and Strauss (2006) for more complete examples.

  # fit the nonparametric model to give in-sample predictions
  # store the parameters to make population inference later
  res <- ecoNP(Y ~ X, data = reg, n.draws = 50, param = TRUE, verbose = TRUE)
  #summarize the results
  x <- summary(res)
  expect_that(length(x), is_equivalent_to(8))
  expect_true(is.null(x$agg.wtable))
  expect_equal(x$agg.table[1,2], 0.04059766, tolerance = accuracy1)
  expect_equal(x$agg.table[2,3], 0.8129786, tolerance = accuracy1)

  # obtain out-of-sample prediction
  out <- predict(res, verbose = TRUE)
  # summarize the results
  x <- summary(out)
  expect_that(length(x), is_equivalent_to(2))
  expect_true("n.draws" %in% names(x))
  expect_equal(x$W.table[1,3], 0.02617743, tolerance = accuracy1)
  expect_equal(x$W.table[2,1], 0.8137116, tolerance = accuracy1)

  # density plots of the out-of-sample predictions
  # par(mfrow=c(2,1))
  # plot(density(out[,1]), main = "W1")
  # plot(density(out[,2]), main = "W2")
})

if (!donotrun) test_that("tests ecoNP on Robinson census data", {
  # load the Robinson's census data
  data(census)

  # fit the parametric model with contextual effects and N 
  # using the default prior specification

  res1 <- ecoNP(Y ~ X, N = N, context = TRUE, param = TRUE, data = census, n.draws = 25, verbose = TRUE)
  # summarize the results
  x <- summary(res1)
  expect_that(length(x), is_equivalent_to(8))
  expect_false(is.null(x$agg.wtable))
  expect_equal(x$agg.table[1,2], 0.009952511, tolerance = accuracy1)
  expect_equal(x$agg.table[2,3], 0.8690776, tolerance = accuracy1)
  expect_equal(x$agg.wtable[1,2], 0.009508983, tolerance = accuracy1)
  expect_equal(x$agg.wtable[2,3], 0.9005222, tolerance = accuracy1)
 
  # out-of sample prediction 
  pres1 <- predict(res1)
  x <- summary(pres1)
  expect_that(length(x), is_equivalent_to(2))
  expect_true("n.draws" %in% names(x))
  expect_equal(x$W.table[1,3], 0.1333375, tolerance = accuracy1)
  expect_equal(x$W.table[2,1], 0.8434944, tolerance = accuracy1)
})