if (FALSE) { if (requireNamespace("tinyVAST", quietly = TRUE)) { library("tinyVAST", warn.conflicts = FALSE) TOL <- 1e-5 get_sdmTMB_pars <- function(x) { c( as.list(x$sd_report, "Estimate"), as.list(x$sd_report, "Estimate", report = TRUE) ) } test_that("tinyVAST/sdmTMB Tweedie spatiotemporal IID models and index area integration match", { skip_on_cran() mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 18) system.time({ fit_sd <- sdmTMB( density ~ 0 + as.factor(year), data = pcod, mesh = mesh, family = tweedie(), time = "year", control = sdmTMBcontrol(multiphase = FALSE) ) }) ps <- get_sdmTMB_pars(fit_sd) system.time({ fit_tv <- tinyVAST( density ~ 0 + factor(year), dsem = "", sem = "", data = pcod, family = tweedie(), time_column = "year", space_columns = c("X", "Y"), spatial_graph = mesh$mesh, control = tinyVASTcontrol(newton_loops = 1) ) }) pt <- as.list(fit_tv$sdrep, "Estimate") expect_equal(pt$alpha_j, ps$b_j, tolerance = TOL) expect_equal(pt$log_sigma[2], ps$thetaf, tolerance = TOL) expect_equal(pt$log_sigma[1], ps$ln_phi, tolerance = TOL) expect_equal(pt$beta_z, ps$sigma_E[1, 1], tolerance = TOL) expect_equal(pt$theta_z, ps$sigma_O[1, 1], tolerance = TOL) expect_equal(pt$log_kappa, ps$ln_kappa[1, 1], tolerance = TOL) g <- replicate_df(qcs_grid, "year", unique(pcod$year)) p <- predict(fit_sd, newdata = g, return_tmb_object = TRUE) system.time({ is <- get_index(p, bias_correct = TRUE) }) # system.time({ # is2 <- lapply(unique(g$year), \(x) { # pp <- predict(fit_sd, newdata = subset(g, year == x), return_tmb_object = TRUE) # get_index(pp, bias_correct = TRUE) # }) # is2 <- do.call(rbind, is2) # }) system.time({ g$var <- "response" it <- lapply(unique(g$year), \(x) integrate_output(fit_tv, newdata = subset(g, year == x), apply.epsilon = TRUE)) it <- do.call(rbind, it) |> as.data.frame() }) expect_equal(it$`Est. (bias.correct)`, is$est, tolerance = TOL) }) fit_sd <- sdmTMB( density ~ 0 + as.factor(year), data = pcod, family = delta_lognormal(type = "poisson-link"), spatial = "off", spatiotemporal = "off" ) ps <- get_pars(fit_sd) suppressWarnings({ fit_tv <- tinyVAST( density ~ 0 + factor(year), delta_options = list(delta_formula = ~ 0 + factor(year)), data = pcod, family = delta_lognormal(type = "poisson-link"), control = tinyVASTcontrol(newton_loops = 1) ) }) pt <- as.list(fit_tv$sdrep, "Estimate") compare_models_delta <- function(mtv, msd) { ps <- get_pars(msd) pt <- as.list(mtv$sdrep, "Estimate") expect_equal(logLik(msd), logLik(mtv), ignore_attr = TRUE) expect_equal(pt$alpha_j, ps$b_j, tolerance = TOL) expect_equal(pt$alpha2_j, ps$b_j2, tolerance = TOL) } test_that("tinyVAST/sdmTMB delta_lognormal(type = 'poisson-link') models match", { skip_on_cran() compare_models_delta(fit_tv, fit_sd) expect_equal(pt$log_sigma, ps$ln_phi[2], tolerance = TOL) }) test_that("tinyVAST/sdmTMB delta_gamma() models match", { skip_on_cran() fit_sd <- update(fit_sd, family = delta_gamma()) fit_tv <- update(fit_tv, family = delta_gamma()) compare_models_delta(fit_tv, fit_sd) ps <- get_pars(fit_sd) pt <- as.list(fit_tv$sdrep, "Estimate") expect_equal(1 / (exp(pt$log_sigma))^2, exp(ps$ln_phi[2]), tolerance = TOL) }) test_that("tinyVAST/sdmTMB delta_gamma(type = 'poisson-link')) models match", { skip_on_cran() fit_sd <- update(fit_sd, family = delta_gamma(type = "poisson-link")) fit_tv <- update(fit_tv, family = delta_gamma(type = "poisson-link")) compare_models_delta(fit_tv, fit_sd) ps <- get_pars(fit_sd) pt <- as.list(fit_tv$sdrep, "Estimate") expect_equal(1 / (exp(pt$log_sigma))^2, exp(ps$ln_phi[2]), tolerance = TOL) }) } }