tests/testthat/test-diag-derive.R

In vlandau/tempo: Temporally-Explicit Models for Phenological Observations

context("tempo derived quantities")
test_that(
  "tests for derived quantites", {
    ### Simulate some data
    library(tempo)
    data_re <- tempo_simulate(n = 20, t = 50,
                              censor_intensity = 3,
                              beta_mu = c(-5, 0.7, -1, 0.3),
                              random_betas = c(1, 3),
                              beta_sd = c(0.23, 0.3),
                              covariates = NULL,
                              n_group = 3,
                              rho = 0)
    
    data_fe <- tempo_simulate(n = 20, t = 50,
                              censor_intensity = 3,
                              beta_mu = c(-5, 0.5, -1, 0.3))
    ### Run sampler for each
    draws_re <- tempo_mcmc(y = data_re$y,
                           covariates = data_re$covariates, 
                           beta_start = c(-5, 0.5, -1, 0.3),
                           n_chains = 1,
                           n_burnin = 50,
                           n_iter = 40,
                           group_ids = data_re$j,
                           sd_eps_start = c(0.05, 0.05),
                           correlated = FALSE,
                           monitor_random_effects = TRUE,
                           random_effects = c("intercept", "covariate1"),
                           pp_checks = "mean")
    
    draws_fe <- tempo_mcmc(y = data_fe$y,
                           covariates = data_fe$covariates, 
                           beta_start = c(-5, 0.5, -1, 0.3),
                           n_chains = 3,
                           n_burnin = 50,
                           n_iter = 40,
                           pp_checks = "variance")
    
    ##### tempo_derive tests #####
    ### calculate derived quantities for both using functions
    dquant_re <- tempo_derive(draws_re, 
                              quantity = c("probabilities",
                                           "pmf", "cdf",
                                           "mean", "sd"),
                              obs_idx = c(1,20),
                              n_cores = 2)
    
    dquant_fe <- tempo_derive(draws_fe, 
                              quantity = c("probabilities",
                                           "pmf", "cdf",
                                           "mean", "sd"),
                              obs_idx = c(11),
                              parallelize = FALSE,
                              n_cores = NULL)
    
    ### calc derived quantities manually to check functions are correct
    ## random effects
    re_cov_array <- do.call(abind, c(data_re$covariates, along = 3))
    # for observation 1 (group 1), iteration 10
    betas_re_1_10 <- c(draws_re$chain_1$mean_beta_intercept[10] +
                        draws_re$chain_1$eps_intercept_1[10],
                      draws_re$chain_1$mean_beta_covariate1[10] +
                        draws_re$chain_1$eps_covariate1_1[10],
                      draws_re$chain_1$beta_covariate2[10],
                      draws_re$chain_1$beta_covariate3[10])
    probs_re_1_10 <- boot::inv.logit(
      betas_re_1_10 %*% rbind(1, t(re_cov_array[1,,]))
    )
  
    pmf_re_1_10 <- dtvgeom(1:51, probs_re_1_10)
    cdf_re_1_10 <- cumsum(pmf_re_1_10)
    mean_re_1_10 <- tvgeom_mean(probs_re_1_10)
    sd_re_1_10 <- sqrt(tvgeom_var(probs_re_1_10))

    expect_equal(unname(dquant_re$probabilities[10, , "obs_1"]),
                     as.vector(probs_re_1_10))
    expect_equal(unname(dquant_re$pmf[10, , "obs_1"]),
                 as.vector(pmf_re_1_10))
    expect_equal(unname(dquant_re$cdf[10, , "obs_1"]),
                 as.vector(cdf_re_1_10))
    expect_equal(dquant_re$mean[10, "obs_1"],
                 mean_re_1_10)
    expect_equal(dquant_re$sd[10, "obs_1"],
                 sd_re_1_10)
    
    # for observation 20 (group 3), iteration 10
    betas_re_3_10 <- c(draws_re$chain_1$mean_beta_intercept[10] +
                         draws_re$chain_1$eps_intercept_3[10],
                       draws_re$chain_1$mean_beta_covariate1[10] +
                         draws_re$chain_1$eps_covariate1_3[10],
                       draws_re$chain_1$beta_covariate2[10],
                       draws_re$chain_1$beta_covariate3[10])
    probs_re_3_10 <- boot::inv.logit(
      betas_re_3_10 %*% rbind(1, t(re_cov_array[20, , ]))
    )
    
    pmf_re_3_10 <- dtvgeom(1:51, probs_re_3_10)
    cdf_re_3_10 <- cumsum(pmf_re_3_10)
    mean_re_3_10 <- tvgeom_mean(probs_re_3_10)
    sd_re_3_10 <- sqrt(tvgeom_var(probs_re_3_10))
    
    expect_equal(unname(dquant_re$probabilities[10, , "obs_20"]),
                 as.vector(probs_re_3_10))
    expect_equal(unname(dquant_re$pmf[10, , "obs_20"]),
                 as.vector(pmf_re_3_10))
    expect_equal(unname(dquant_re$cdf[10, , "obs_20"]),
                 as.vector(cdf_re_3_10))
    expect_equal(dquant_re$mean[10, "obs_20"],
                 mean_re_3_10)
    expect_equal(dquant_re$sd[10, "obs_20"],
                 sd_re_3_10)
    
    ## fixed effects
    fe_cov_array <- do.call(abind, c(data_fe$covariates, along = 3))
    # for observation 11, iteration 10
    betas_fe_11_10 <- c(draws_fe$chain_1$beta_intercept[10],
                        draws_fe$chain_1$beta_covariate1[10],
                        draws_fe$chain_1$beta_covariate2[10],
                        draws_fe$chain_1$beta_covariate3[10])
    probs_fe_11_10 <- boot::inv.logit(
      betas_fe_11_10 %*% rbind(1, t(fe_cov_array[11,,]))
    )
    
    pmf_fe_11_10 <- dtvgeom(1:51, probs_fe_11_10)
    cdf_fe_11_10 <- cumsum(pmf_fe_11_10)
    mean_fe_11_10 <- tvgeom_mean(probs_fe_11_10)
    sd_fe_11_10 <- sqrt(tvgeom_var(probs_fe_11_10))
    
    expect_equal(unname(dquant_fe$probabilities[10, , "obs_11"]),
                 as.vector(probs_fe_11_10))
    expect_equal(unname(dquant_fe$pmf[10, , "obs_11"]),
                 as.vector(pmf_fe_11_10))
    expect_equal(unname(dquant_fe$cdf[10, , "obs_11"]),
                 as.vector(cdf_fe_11_10))
    expect_equal(dquant_fe$mean[10, "obs_11"],
                 mean_fe_11_10)
    expect_equal(dquant_fe$sd[10, "obs_11"],
                 sd_fe_11_10)
    
    
    ## Test error throws
    expect_error(tempo_derive(draws_fe, quantity = "abc", obs_idx = 3))
    
    ##### tempo_diagnostics tests #####
    
    ## Syntax checks
    expect_error({tempo_dic(draws_re, n_cores = 2)}, NA)
    expect_error({tempo_dic(draws_fe, parallelize = FALSE)}, NA)
    expect_error(tempo_psrf(draws_re))
    expect_error({tempo_psrf(draws_fe)}, NA)
    
  }
)