tests/testthat/test_calc_stats.R

In kenkellner/jagsUI2: A Wrapper Around 'rjags' to Streamline 'JAGS' Analyses

context("Test calculation accuracy")

test_that("Summary calculations match reference",{
  skip_on_cran()
  data(longley)
  jags_data <- list(gnp=longley$GNP, employed=longley$Employed, 
             n=length(longley$Employed))
  model_file <- tempfile()
  writeLines("
  model{
  #Likelihood
  for (i in 1:n){ 
  employed[i] ~ dnorm(mu[i], tau)     
  mu[i] <- alpha + beta*gnp[i]
  }
  #Priors
  alpha ~ dnorm(0, 0.00001)
  beta ~ dnorm(0, 0.00001)
  sigma ~ dunif(0,1000)
  tau <- pow(sigma,-2)
  delta ~ dnorm(-1,0.01) # to test overlap0
  }", con=model_file)
  params <- c('alpha','beta','sigma','mu','delta')
  n_chains <- 3; n_iter <<- 1000; n_warmup <<- 500
  set.seed(123)
  inits <- vector(length=3,"list")
  for (i in 1:n_chains){
    inits[[i]]$.RNG.name<-"base::Mersenne-Twister"
    inits[[i]]$.RNG.seed<-floor(stats::runif(1,1,100000))
  }
  example_summary <- readRDS(file.path("example_summary.rds"))
  out <- jagsUI(model_file,jags_data,inits=inits,params,
              n_chains=n_chains,n_iter=n_iter,n_warmup=n_warmup, quiet=T)
  expect_equal(out$summary, example_summary)
})

test_that("calc_stats handles NAs correctly",{
  samples <- readRDS('coda_samples.rds')
  samples[[1]][1,'sigma'] <- NA
  stats_1na <- jagsUI2:::calc_stats(samples)
  expect_true(all(is.na(stats_1na['sigma',])))
  expect_false(any(is.na(stats_1na['alpha',])))

  samples[[1]][,'sigma'] <- NA
  stats_allna <- jagsUI2:::calc_stats(samples)
  expect_true(all(is.na(stats_allna['sigma',])))
  expect_false(any(is.na(stats_allna['alpha',])))

})

test_that("Inf in samples handled correctly",{
  samples <- readRDS('coda_samples.rds')
  samples[[1]][1,'sigma'] <- Inf
  stats_inf <- jagsUI2:::calc_stats(samples)
  expect_equal(stats_inf['sigma',1],Inf)
  expect_true(all(is.na(stats_inf['sigma',2:9])))
  expect_false(any(is.na(stats_inf['alpha',])))
})

test_that("Calculation of pD/DIC works", {
  samples <- readRDS('coda_samples.rds')
  expect_equal(calc_DIC(samples), c(pD=6.660906,DIC=40.712014), tol=1e-4)
  no_dev <- remove_params(samples, 'deviance')
  expect_true(all(is.na(calc_DIC(no_dev))))
  samp_na <- samples
  ind <- which_params('deviance',param_names(samples))
  samp_na[[1]][1,ind] <- NA 
  expect_true(all(is.na(calc_DIC(samp_na))))
  samp_inf <- samples
  samp_inf[[1]][1,ind] <- Inf
  expect_true(all(is.na(calc_DIC(samp_na))))
  samp_inf[[1]][1,ind] <- -Inf
  expect_true(all(is.na(calc_DIC(samp_na))))
})