tests/testthat/test-2-fit-less-basic.R
In pbs-assess/sdmTMB: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
test_that("Regularization works", {
skip_on_cran()
local_edition(2)
d <- subset(pcod, year >= 2015)
d$depth_scaled <- as.numeric(scale(d$depth_scaled))
m1 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"))
# Bypassing via NAs:
m2 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"),
priors = sdmTMBpriors(b = normal(c(NA, NA, NA), c(NA, NA, NA))))
expect_equal(m1$sd_report, m2$sd_report)
# Ridge regression on depth term:
m2 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"),
priors = sdmTMBpriors(b = normal(c(0, NA, NA), c(1, NA, NA))))
b1 <- tidy(m1)
b2 <- tidy(m2)
expect_lt(b2$estimate[1], b1$estimate[2])
expect_lt(b2$std.error[1], b1$std.error[2])
})
test_that("A time-varying model fits and predicts appropriately", {
skip_on_cran()
local_edition(2)
# SEED <- 42
# set.seed(SEED)
# x <- stats::runif(60, -1, 1)
# y <- stats::runif(60, -1, 1)
# initial_betas <- 0.5
# range <- 0.5
# sigma_O <- 0
# sigma_E <- 0.1
# phi <- 0.1
# sigma_V <- 0.3
# loc <- data.frame(x = x, y = y)
# spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)
#
# s <- sdmTMB_sim(
# x = x, y = y, mesh = spde, range = range,
# betas = initial_betas, time_steps = 12L, sigma_V = sigma_V,
# phi = phi, sigma_O = sigma_O, sigma_E = sigma_E,
# seed = SEED
# )
spde <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
m <- sdmTMB(data = pcod, formula = density ~ 0, spatial = TRUE,
time_varying = ~ 1, time = "year", mesh = spde, family = tweedie(),
spatiotemporal = "off")
expect_equal(exp(m$model$par["ln_tau_V"])[[1]], 0.5971512, tolerance = 0.001)
tidy(m, effects = "ran_par")
# b_t <- dplyr::group_by(s, time) %>%
# dplyr::summarize(b_t = unique(b), .groups = "drop") %>%
# dplyr::pull(b_t)
# r <- m$tmb_obj$report()
# b_t_fit <- r$b_rw_t[,,1]
# plot(b_t, b_t_fit, asp = 1);abline(a = 0, b = 1)
# expect_equal(mean((b_t- b_t_fit)^2), 0, tolerance = 1e-4)
p <- predict(m, newdata = NULL)
# plot(p$est, s$observed, asp = 1);abline(a = 0, b = 1)
# expect_equal(mean((p$est - s$observed)^2), 0, tolerance = 0.01)
cols <- c("est", "est_non_rf", "est_rf", "omega_s")
p_nd <- predict(m, newdata = pcod)
expect_equal(p[,cols], p_nd[,cols], tolerance = 1e-4)
})
test_that("Year indexes get created correctly", {
expect_identical(make_year_i(c(1, 2, 3)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1L, 2L, 3L)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1L, 2L, 4L)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1, 2, 4, 2)), c(0L, 1L, 2L, 1L))
expect_identical(make_year_i(c(1, 4, 2)), c(0L, 2L, 1L))
})
test_that("A spatially varying coefficient model fits", {
skip_on_cran()
d <- subset(pcod, year >= 2011) # subset for speed
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30)
d$scaled_year <- (d$year - mean(d$year)) / sd(d$year)
m <- sdmTMB(density ~ depth_scaled, data = d,
mesh = pcod_spde, family = tweedie(link = "log"),
spatial_varying = ~ 0 + scaled_year, time = "year")
expect_true(all(!is.na(summary(m$sd_report)[,"Std. Error"])))
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$scaled_year <- scale(nd$year)
p <- predict(m, newdata = subset(nd, year >= 2011))
})
test_that("SPDE as generated by make_mesh is consistent", {
set.seed(42)
skip_on_cran()
set.seed(1)
local_edition(2)
x <- stats::runif(70, -1, 1)
y <- stats::runif(70, -1, 1)
loc <- data.frame(x = x, y = y)
spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)
# spot check first 10 locations
loc_xy <- matrix(c(
-0.468982674,-0.321854124,
-0.255752201,0.678880700,
0.145706727,-0.306633022,
0.816415580,-0.332450138,
-0.596636138,-0.047297510,
0.796779370,0.784396672,
0.889350537,0.728678941,
0.321595585,-0.220020913,
0.258228088,0.554641398,
-0.876427459,0.921235994), ncol = 2, byrow = TRUE)
expect_equal(c(spde$loc_xy[1:10,]), c(loc_xy[1:10,]), tolerance = 1e-8)
# Depends on stable or testing INLA!?
# test the n
# expect_equal(spde$spde$param.inla$n, 141.0, tolerance = 1e-8)
# expect_equal(spde$spde$param.inla$n, 150)
# test M0 values - first 10
# expect_equal(c(0.06889138, 0.04169878, 0.01547313, 0.05734409, 0.03271476, 0.03672571, 0.05200510, 0.01293629, 0.03631651, 0.04446844),
# as.numeric(spde$spde$param.inla$M0)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
# # test M1 values - first 10
# expect_equal(c(3.602458, 5.105182, 4.571970, 4.353246, 4.496775, 4.802626, 4.505773, 4.935963, 4.419817, 5.169096),
# as.numeric(spde$spde$param.inla$M1)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
# # test M2 values - first 10
# expect_equal(c(253.1738,814.3858,1754.2543,401.3752,865.0118,754.9644,543.8968,2425.8404,729.7511,1014.0049),
# as.numeric(spde$spde$param.inla$M2)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
})
# test_that("The `map` argument works.", {
# skip_on_cran()
# skip_on_ci()
# d <- subset(pcod, year >= 2013)
# pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
#
# m <- sdmTMB(density ~ 0 + as.factor(year),
# data = d, time = "year", mesh = pcod_spde,
# family = tweedie(link = "log")
# )
# p <- predict(m)
#
# .map <- m$tmb_map
# .map$ln_tau_O <- as.factor(NA)
# .map$ln_tau_E <- as.factor(NA)
# .map$omega_s <- factor(rep(NA, length(m$tmb_params$omega_s)))
# .map$epsilon_st <- factor(rep(NA, length(m$tmb_params$epsilon_st)))
# .map$ln_kappa <- as.factor(NA)
#
# m.map <- sdmTMB(density ~ 0 + as.factor(year),
# data = d, time = "year", mesh = pcod_spde,
# family = tweedie(link = "log"), map = .map
# )
# p.map <- predict(m.map)
#
# expect_true(all(p.map$omega_s == 0))
# expect_true(all(p.map$epsilon_st == 0))
# expect_true(!identical(p$est, p.map$est))
# })
test_that("profile option works", {
skip_on_cran()
d <- subset(pcod, year == 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 10)
mp <- sdmTMB(density ~ depth_scaled + depth_scaled2,
data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = TRUE),
family = tweedie(link = "log")
)
m <- sdmTMB(density ~ depth_scaled + depth_scaled2,
data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = FALSE),
family = tweedie(link = "log")
)
expect_lt(length(mp$model$par), length(m$model$par))
bp <- tidy(mp)
b <- tidy(m)
expect_equal(bp$estimate, b$estimate, tolerance = 0.05)
})
test_that("The mapping off spatial and spatiotemporal fields works.", {
skip_on_cran()
d <- subset(pcod, year >= 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
m <- sdmTMB(density ~ 0 + as.factor(year),
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log")
)
p <- predict(m)
m.map <- sdmTMB(density ~ 0 + as.factor(year),
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log"), spatial = "off", spatiotemporal = "off"
)
p.map <- predict(m.map)
expect_true(is.na(m.map$tmb_map$ln_tau_E))
expect_true(is.na(m.map$tmb_map$ln_tau_O))
expect_true(!identical(p$est, p.map$est))
expect_true(length(unique(p.map$est)) == length(unique(d$year)))
# basic lm test:
dpos <- d[d$density > 0, ]
pcod_spde <- make_mesh(dpos, c("X", "Y"), cutoff = 50)
m.sdmTMB.map <- sdmTMB(log(density) ~ depth, data = dpos,
family = gaussian(), spatial = "off", spatiotemporal = "off",
mesh = pcod_spde)
m.stats.glm <- stats::lm(log(density) ~ depth, data = dpos)
m.glmmTMB <- glmmTMB::glmmTMB(log(density) ~ depth, data = dpos)
.t <- tidy(m.sdmTMB.map)
expect_equal(.t$estimate, as.numeric(coef(m.stats.glm)), tolerance = 1e-5)
expect_equal(exp(m.sdmTMB.map$model$par[["ln_phi"]]),
glmmTMB::sigma(m.glmmTMB), tolerance = 1e-5)
# Bernoulli:
pcod_binom <- pcod_2011
pcod_binom$present <- ifelse(pcod_binom$density > 0, 1L, 0L)
pcod_spde <- make_mesh(pcod_binom, c("X", "Y"), cutoff = 50)
m.sdmTMB.map <- sdmTMB(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"), spatial = "off", spatiotemporal = "off", mesh = pcod_spde)
m.stats.glm <- stats::glm(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"))
m.glmmTMB <- glmmTMB::glmmTMB(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"))
b1 <- as.numeric(unclass(glmmTMB::fixef(m.glmmTMB))$cond)
b2 <- tidy(m.sdmTMB.map)$estimate
b3 <- as.numeric(coef(m.stats.glm))
expect_equal(b1, b2, tolerance = 1e-6)
expect_equal(b3, b2, tolerance = 1e-6)
})
test_that("Large coordinates cause a warning.", {
skip_on_cran()
d <- subset(pcod, year == 2017)
d$X <- d$X * 1000
d$Y <- d$Y * 1000
expect_warning(pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 70))
})
test_that("Random walk fields work", {
skip_on_cran()
d <- subset(pcod, year >= 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 15)
m_rw <- sdmTMB(density ~ 1,
data = d, time = "year", mesh = pcod_spde,
control = sdmTMBcontrol(newton_loops = 1),
family = tweedie(link = "log"), spatiotemporal = "RW"
)
expect_identical(m_rw$tmb_data$rw_fields, TRUE)
p_rw <- predict(m_rw)
# close to AR1 with high rho:
m_ar1 <- sdmTMB(density ~ 1,
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log"), spatiotemporal = "AR1",
control = sdmTMBcontrol(
newton_loops = 1,
start = list(ar1_phi = qlogis((0.98 + 1) / 2)),
map = list(ar1_phi = factor(NA)))
)
expect_identical(m_ar1$tmb_data$ar1_fields, TRUE)
expect_identical(m_ar1$tmb_data$rw_fields, FALSE)
p_ar1 <- predict(m_ar1)
expect_gt(stats::cor(p_ar1$est, p_rw$est), 0.9)
})
test_that("start works", {
skip_on_cran()
expect_error({
m2 <- sdmTMB(density ~ poly(depth, 2),
data = pcod_2011,
mesh = pcod_mesh_2011, family = tweedie(),
control = sdmTMBcontrol(start = list(ln_kappa = -1.78)))
}, regexp = "kappa")
expect_message({
m2 <- sdmTMB(density ~ poly(depth, 2),
data = pcod_2011,
mesh = pcod_mesh_2011, family = tweedie(),
control = sdmTMBcontrol(start = list(ln_kappa = matrix(c(-1.78, -1.78), ncol = 1L))))
}, regexp = "ln_kappa")
})
test_that("Multiple SVC works", {
skip_on_cran()
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)
pcod$syear <- as.numeric(scale(pcod$year))
fit <- sdmTMB(
density ~ 1,
spatial_varying = ~ 0 + syear + depth_scaled,
mesh = mesh,
family = delta_gamma(),
data = pcod
)
fit$sd_report
s <- as.list(fit$sd_report, "Std. Error")
expect_true(sum(is.na(s$b_j)) == 0L)
fit
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$syear <- as.numeric(scale(nd$year))
p <- predict(fit, newdata = nd)
# p <- predict(fit, newdata = NULL)
})
test_that("More esoteric prediction options work", {
skip_on_cran()
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = delta_gamma()
)
expect_error(p <- predict(fit, nsim = 5, model = "aaa"),
regexp = "model"
)
expect_error(p <- predict(fit, nsim = 5, model = 99),
regexp = "model"
)
p <- predict(fit, nsim = 5, model = NA)
head(p)
p <- predict(fit, nsim = 5, model = 1)
head(p)
# p <- predict(fit, nsim = 5, model = 1, type = "response")
# head(p)
# expect_true(all(p <= 1 & p >= 0))
p <- predict(fit, nsim = 5, model = 2)
head(p)
# p <- predict(fit, nsim = 5, model = 2, type = "response")
# head(p)
# expect_true(all(p > 0))
p <- predict(fit,
nsim = 5, model = 2,
return_tmb_report = TRUE
)
expect_length(p, 5L)
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie()
)
p <- predict(fit, nsim = 5)
# p <- predict(fit, nsim = 5, type = "response")
# head(p)
# expect_true(all(p > 0))
})
test_that("update() works", {
skip_on_cran()
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie()
)
fit2 <- update(fit)
expect_equal(fit$model, fit2$model)
})
test_that("Irregular time gets detected", {
skip_on_cran()
mesh <- make_mesh(pcod, c("X", "Y"), n_knots = 50)
expect_message({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
do_fit = FALSE
)
}, regexp = "irregular time")
expect_message({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
do_fit = FALSE
)
}, regexp = "c\\(2006, 2008, 2010, 2012, 2014, 2016\\)")
expect_silent({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
do_fit = FALSE
)
})
})
test_that("find_missing_time works", {
expect_identical(find_missing_time(c(1, 3, 4)), 2)
expect_identical(find_missing_time(c(1, 3, 4, 5, 7)), c(2, 6))
expect_identical(find_missing_time(c(1, 2, 3)), numeric(0))
expect_identical(find_missing_time(c(1, 4, 3)), 2)
expect_identical(find_missing_time(c(1L, 4L, 3L)), 2L)
})
test_that("offset() throws an error", {
expect_error({
fit <- sdmTMB(
density ~ 1 + offset(depth),
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie(link = "log")
)
}, regexp = "offset")
expect_error({
fit <- sdmTMB(
density ~ 1 + offset(log(depth)),
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie(link = "log")
)
}, regexp = "offset")
})
test_that("sdmTMB throws error on AR1/RW with non-numeric time", {
skip_on_cran()
pcod_2011$fyear <- as.factor(pcod_2011$year)
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "iid", data = pcod_2011,
mesh = pcod_mesh_2011)
expect_true(inherits(fit, "sdmTMB"))
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "ar1", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "rw", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
pcod_2011$cyear <- as.character(pcod_2011$year)
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "iid", data = pcod_2011,
mesh = pcod_mesh_2011)
expect_true(inherits(fit, "sdmTMB"))
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "ar1", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "rw", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
fit <- sdmTMB(
present ~ 0 + depth,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "year",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
)
expect_error(
fit <- sdmTMB(
present ~ 0,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "fyear",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
))
expect_error(
fit <- sdmTMB(
present ~ 0,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "cyear",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
))
})
test_that("Time with an NA gets flagged", {
skip_on_cran()
d <- pcod_2011
d$yr <- as.character(d$year)
d$yr[2] <- NA
expect_error({
m <- sdmTMB(density ~ 1, time = "yr", spatial = "off", spatiotemporal = "off", data = d)
}, regexp = "time")
})
test_that("Prediction outside fitted coordinates gets warned about #285", {
fit <- sdmTMB(
present ~ 1,
family = binomial(),
data = pcod_2011,
mesh = pcod_mesh_2011
)
nd <- qcs_grid
range(nd$X)
nd$X <- nd$X * 10
range(nd$X)
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$X <- nd$X / 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$Y <- nd$Y / 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$Y <- nd$Y * 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
})
pbs-assess/sdmTMB documentation built on May 17, 2024, 11:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
test_that("Regularization works", {
skip_on_cran()
local_edition(2)
d <- subset(pcod, year >= 2015)
d$depth_scaled <- as.numeric(scale(d$depth_scaled))
m1 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"))
# Bypassing via NAs:
m2 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"),
priors = sdmTMBpriors(b = normal(c(NA, NA, NA), c(NA, NA, NA))))
expect_equal(m1$sd_report, m2$sd_report)
# Ridge regression on depth term:
m2 <- sdmTMB(data = d,
formula = density ~ 0 + depth_scaled + as.factor(year),
mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
family = tweedie(link = "log"),
priors = sdmTMBpriors(b = normal(c(0, NA, NA), c(1, NA, NA))))
b1 <- tidy(m1)
b2 <- tidy(m2)
expect_lt(b2$estimate[1], b1$estimate[2])
expect_lt(b2$std.error[1], b1$std.error[2])
})
test_that("A time-varying model fits and predicts appropriately", {
skip_on_cran()
local_edition(2)
# SEED <- 42
# set.seed(SEED)
# x <- stats::runif(60, -1, 1)
# y <- stats::runif(60, -1, 1)
# initial_betas <- 0.5
# range <- 0.5
# sigma_O <- 0
# sigma_E <- 0.1
# phi <- 0.1
# sigma_V <- 0.3
# loc <- data.frame(x = x, y = y)
# spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)
#
# s <- sdmTMB_sim(
# x = x, y = y, mesh = spde, range = range,
# betas = initial_betas, time_steps = 12L, sigma_V = sigma_V,
# phi = phi, sigma_O = sigma_O, sigma_E = sigma_E,
# seed = SEED
# )
spde <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
m <- sdmTMB(data = pcod, formula = density ~ 0, spatial = TRUE,
time_varying = ~ 1, time = "year", mesh = spde, family = tweedie(),
spatiotemporal = "off")
expect_equal(exp(m$model$par["ln_tau_V"])[[1]], 0.5971512, tolerance = 0.001)
tidy(m, effects = "ran_par")
# b_t <- dplyr::group_by(s, time) %>%
# dplyr::summarize(b_t = unique(b), .groups = "drop") %>%
# dplyr::pull(b_t)
# r <- m$tmb_obj$report()
# b_t_fit <- r$b_rw_t[,,1]
# plot(b_t, b_t_fit, asp = 1);abline(a = 0, b = 1)
# expect_equal(mean((b_t- b_t_fit)^2), 0, tolerance = 1e-4)
p <- predict(m, newdata = NULL)
# plot(p$est, s$observed, asp = 1);abline(a = 0, b = 1)
# expect_equal(mean((p$est - s$observed)^2), 0, tolerance = 0.01)
cols <- c("est", "est_non_rf", "est_rf", "omega_s")
p_nd <- predict(m, newdata = pcod)
expect_equal(p[,cols], p_nd[,cols], tolerance = 1e-4)
})
test_that("Year indexes get created correctly", {
expect_identical(make_year_i(c(1, 2, 3)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1L, 2L, 3L)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1L, 2L, 4L)), c(0L, 1L, 2L))
expect_identical(make_year_i(c(1, 2, 4, 2)), c(0L, 1L, 2L, 1L))
expect_identical(make_year_i(c(1, 4, 2)), c(0L, 2L, 1L))
})
test_that("A spatially varying coefficient model fits", {
skip_on_cran()
d <- subset(pcod, year >= 2011) # subset for speed
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30)
d$scaled_year <- (d$year - mean(d$year)) / sd(d$year)
m <- sdmTMB(density ~ depth_scaled, data = d,
mesh = pcod_spde, family = tweedie(link = "log"),
spatial_varying = ~ 0 + scaled_year, time = "year")
expect_true(all(!is.na(summary(m$sd_report)[,"Std. Error"])))
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$scaled_year <- scale(nd$year)
p <- predict(m, newdata = subset(nd, year >= 2011))
})
test_that("SPDE as generated by make_mesh is consistent", {
set.seed(42)
skip_on_cran()
set.seed(1)
local_edition(2)
x <- stats::runif(70, -1, 1)
y <- stats::runif(70, -1, 1)
loc <- data.frame(x = x, y = y)
spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)
# spot check first 10 locations
loc_xy <- matrix(c(
-0.468982674,-0.321854124,
-0.255752201,0.678880700,
0.145706727,-0.306633022,
0.816415580,-0.332450138,
-0.596636138,-0.047297510,
0.796779370,0.784396672,
0.889350537,0.728678941,
0.321595585,-0.220020913,
0.258228088,0.554641398,
-0.876427459,0.921235994), ncol = 2, byrow = TRUE)
expect_equal(c(spde$loc_xy[1:10,]), c(loc_xy[1:10,]), tolerance = 1e-8)
# Depends on stable or testing INLA!?
# test the n
# expect_equal(spde$spde$param.inla$n, 141.0, tolerance = 1e-8)
# expect_equal(spde$spde$param.inla$n, 150)
# test M0 values - first 10
# expect_equal(c(0.06889138, 0.04169878, 0.01547313, 0.05734409, 0.03271476, 0.03672571, 0.05200510, 0.01293629, 0.03631651, 0.04446844),
# as.numeric(spde$spde$param.inla$M0)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
# # test M1 values - first 10
# expect_equal(c(3.602458, 5.105182, 4.571970, 4.353246, 4.496775, 4.802626, 4.505773, 4.935963, 4.419817, 5.169096),
# as.numeric(spde$spde$param.inla$M1)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
# # test M2 values - first 10
# expect_equal(c(253.1738,814.3858,1754.2543,401.3752,865.0118,754.9644,543.8968,2425.8404,729.7511,1014.0049),
# as.numeric(spde$spde$param.inla$M2)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
})
# test_that("The `map` argument works.", {
# skip_on_cran()
# skip_on_ci()
# d <- subset(pcod, year >= 2013)
# pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
#
# m <- sdmTMB(density ~ 0 + as.factor(year),
# data = d, time = "year", mesh = pcod_spde,
# family = tweedie(link = "log")
# )
# p <- predict(m)
#
# .map <- m$tmb_map
# .map$ln_tau_O <- as.factor(NA)
# .map$ln_tau_E <- as.factor(NA)
# .map$omega_s <- factor(rep(NA, length(m$tmb_params$omega_s)))
# .map$epsilon_st <- factor(rep(NA, length(m$tmb_params$epsilon_st)))
# .map$ln_kappa <- as.factor(NA)
#
# m.map <- sdmTMB(density ~ 0 + as.factor(year),
# data = d, time = "year", mesh = pcod_spde,
# family = tweedie(link = "log"), map = .map
# )
# p.map <- predict(m.map)
#
# expect_true(all(p.map$omega_s == 0))
# expect_true(all(p.map$epsilon_st == 0))
# expect_true(!identical(p$est, p.map$est))
# })
test_that("profile option works", {
skip_on_cran()
d <- subset(pcod, year == 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 10)
mp <- sdmTMB(density ~ depth_scaled + depth_scaled2,
data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = TRUE),
family = tweedie(link = "log")
)
m <- sdmTMB(density ~ depth_scaled + depth_scaled2,
data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = FALSE),
family = tweedie(link = "log")
)
expect_lt(length(mp$model$par), length(m$model$par))
bp <- tidy(mp)
b <- tidy(m)
expect_equal(bp$estimate, b$estimate, tolerance = 0.05)
})
test_that("The mapping off spatial and spatiotemporal fields works.", {
skip_on_cran()
d <- subset(pcod, year >= 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
m <- sdmTMB(density ~ 0 + as.factor(year),
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log")
)
p <- predict(m)
m.map <- sdmTMB(density ~ 0 + as.factor(year),
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log"), spatial = "off", spatiotemporal = "off"
)
p.map <- predict(m.map)
expect_true(is.na(m.map$tmb_map$ln_tau_E))
expect_true(is.na(m.map$tmb_map$ln_tau_O))
expect_true(!identical(p$est, p.map$est))
expect_true(length(unique(p.map$est)) == length(unique(d$year)))
# basic lm test:
dpos <- d[d$density > 0, ]
pcod_spde <- make_mesh(dpos, c("X", "Y"), cutoff = 50)
m.sdmTMB.map <- sdmTMB(log(density) ~ depth, data = dpos,
family = gaussian(), spatial = "off", spatiotemporal = "off",
mesh = pcod_spde)
m.stats.glm <- stats::lm(log(density) ~ depth, data = dpos)
m.glmmTMB <- glmmTMB::glmmTMB(log(density) ~ depth, data = dpos)
.t <- tidy(m.sdmTMB.map)
expect_equal(.t$estimate, as.numeric(coef(m.stats.glm)), tolerance = 1e-5)
expect_equal(exp(m.sdmTMB.map$model$par[["ln_phi"]]),
glmmTMB::sigma(m.glmmTMB), tolerance = 1e-5)
# Bernoulli:
pcod_binom <- pcod_2011
pcod_binom$present <- ifelse(pcod_binom$density > 0, 1L, 0L)
pcod_spde <- make_mesh(pcod_binom, c("X", "Y"), cutoff = 50)
m.sdmTMB.map <- sdmTMB(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"), spatial = "off", spatiotemporal = "off", mesh = pcod_spde)
m.stats.glm <- stats::glm(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"))
m.glmmTMB <- glmmTMB::glmmTMB(present ~ depth, data = pcod_binom,
family = binomial(link = "logit"))
b1 <- as.numeric(unclass(glmmTMB::fixef(m.glmmTMB))$cond)
b2 <- tidy(m.sdmTMB.map)$estimate
b3 <- as.numeric(coef(m.stats.glm))
expect_equal(b1, b2, tolerance = 1e-6)
expect_equal(b3, b2, tolerance = 1e-6)
})
test_that("Large coordinates cause a warning.", {
skip_on_cran()
d <- subset(pcod, year == 2017)
d$X <- d$X * 1000
d$Y <- d$Y * 1000
expect_warning(pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 70))
})
test_that("Random walk fields work", {
skip_on_cran()
d <- subset(pcod, year >= 2013)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 15)
m_rw <- sdmTMB(density ~ 1,
data = d, time = "year", mesh = pcod_spde,
control = sdmTMBcontrol(newton_loops = 1),
family = tweedie(link = "log"), spatiotemporal = "RW"
)
expect_identical(m_rw$tmb_data$rw_fields, TRUE)
p_rw <- predict(m_rw)
# close to AR1 with high rho:
m_ar1 <- sdmTMB(density ~ 1,
data = d, time = "year", mesh = pcod_spde,
family = tweedie(link = "log"), spatiotemporal = "AR1",
control = sdmTMBcontrol(
newton_loops = 1,
start = list(ar1_phi = qlogis((0.98 + 1) / 2)),
map = list(ar1_phi = factor(NA)))
)
expect_identical(m_ar1$tmb_data$ar1_fields, TRUE)
expect_identical(m_ar1$tmb_data$rw_fields, FALSE)
p_ar1 <- predict(m_ar1)
expect_gt(stats::cor(p_ar1$est, p_rw$est), 0.9)
})
test_that("start works", {
skip_on_cran()
expect_error({
m2 <- sdmTMB(density ~ poly(depth, 2),
data = pcod_2011,
mesh = pcod_mesh_2011, family = tweedie(),
control = sdmTMBcontrol(start = list(ln_kappa = -1.78)))
}, regexp = "kappa")
expect_message({
m2 <- sdmTMB(density ~ poly(depth, 2),
data = pcod_2011,
mesh = pcod_mesh_2011, family = tweedie(),
control = sdmTMBcontrol(start = list(ln_kappa = matrix(c(-1.78, -1.78), ncol = 1L))))
}, regexp = "ln_kappa")
})
test_that("Multiple SVC works", {
skip_on_cran()
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)
pcod$syear <- as.numeric(scale(pcod$year))
fit <- sdmTMB(
density ~ 1,
spatial_varying = ~ 0 + syear + depth_scaled,
mesh = mesh,
family = delta_gamma(),
data = pcod
)
fit$sd_report
s <- as.list(fit$sd_report, "Std. Error")
expect_true(sum(is.na(s$b_j)) == 0L)
fit
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$syear <- as.numeric(scale(nd$year))
p <- predict(fit, newdata = nd)
# p <- predict(fit, newdata = NULL)
})
test_that("More esoteric prediction options work", {
skip_on_cran()
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = delta_gamma()
)
expect_error(p <- predict(fit, nsim = 5, model = "aaa"),
regexp = "model"
)
expect_error(p <- predict(fit, nsim = 5, model = 99),
regexp = "model"
)
p <- predict(fit, nsim = 5, model = NA)
head(p)
p <- predict(fit, nsim = 5, model = 1)
head(p)
# p <- predict(fit, nsim = 5, model = 1, type = "response")
# head(p)
# expect_true(all(p <= 1 & p >= 0))
p <- predict(fit, nsim = 5, model = 2)
head(p)
# p <- predict(fit, nsim = 5, model = 2, type = "response")
# head(p)
# expect_true(all(p > 0))
p <- predict(fit,
nsim = 5, model = 2,
return_tmb_report = TRUE
)
expect_length(p, 5L)
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie()
)
p <- predict(fit, nsim = 5)
# p <- predict(fit, nsim = 5, type = "response")
# head(p)
# expect_true(all(p > 0))
})
test_that("update() works", {
skip_on_cran()
fit <- sdmTMB(
density ~ depth_scaled,
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie()
)
fit2 <- update(fit)
expect_equal(fit$model, fit2$model)
})
test_that("Irregular time gets detected", {
skip_on_cran()
mesh <- make_mesh(pcod, c("X", "Y"), n_knots = 50)
expect_message({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
do_fit = FALSE
)
}, regexp = "irregular time")
expect_message({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
do_fit = FALSE
)
}, regexp = "c\\(2006, 2008, 2010, 2012, 2014, 2016\\)")
expect_silent({
fit <- sdmTMB(density ~ 1, time = "year",
data = pcod, mesh = mesh,
family = tweedie(), spatiotemporal = "ar1",
extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
do_fit = FALSE
)
})
})
test_that("find_missing_time works", {
expect_identical(find_missing_time(c(1, 3, 4)), 2)
expect_identical(find_missing_time(c(1, 3, 4, 5, 7)), c(2, 6))
expect_identical(find_missing_time(c(1, 2, 3)), numeric(0))
expect_identical(find_missing_time(c(1, 4, 3)), 2)
expect_identical(find_missing_time(c(1L, 4L, 3L)), 2L)
})
test_that("offset() throws an error", {
expect_error({
fit <- sdmTMB(
density ~ 1 + offset(depth),
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie(link = "log")
)
}, regexp = "offset")
expect_error({
fit <- sdmTMB(
density ~ 1 + offset(log(depth)),
data = pcod_2011, mesh = pcod_mesh_2011,
family = tweedie(link = "log")
)
}, regexp = "offset")
})
test_that("sdmTMB throws error on AR1/RW with non-numeric time", {
skip_on_cran()
pcod_2011$fyear <- as.factor(pcod_2011$year)
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "iid", data = pcod_2011,
mesh = pcod_mesh_2011)
expect_true(inherits(fit, "sdmTMB"))
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "ar1", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "fyear", spatiotemporal = "rw", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
pcod_2011$cyear <- as.character(pcod_2011$year)
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "iid", data = pcod_2011,
mesh = pcod_mesh_2011)
expect_true(inherits(fit, "sdmTMB"))
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "ar1", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
expect_error(
fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
time = "cyear", spatiotemporal = "rw", data = pcod_2011,
mesh = pcod_mesh_2011), regexp = "Time"
)
fit <- sdmTMB(
present ~ 0 + depth,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "year",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
)
expect_error(
fit <- sdmTMB(
present ~ 0,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "fyear",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
))
expect_error(
fit <- sdmTMB(
present ~ 0,
family = binomial(),
spatial = 'off',
time_varying = ~ 1,
time = "cyear",
data = pcod_2011,
spatiotemporal = "off",
mesh = pcod_mesh_2011
))
})
test_that("Time with an NA gets flagged", {
skip_on_cran()
d <- pcod_2011
d$yr <- as.character(d$year)
d$yr[2] <- NA
expect_error({
m <- sdmTMB(density ~ 1, time = "yr", spatial = "off", spatiotemporal = "off", data = d)
}, regexp = "time")
})
test_that("Prediction outside fitted coordinates gets warned about #285", {
fit <- sdmTMB(
present ~ 1,
family = binomial(),
data = pcod_2011,
mesh = pcod_mesh_2011
)
nd <- qcs_grid
range(nd$X)
nd$X <- nd$X * 10
range(nd$X)
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$X <- nd$X / 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$Y <- nd$Y / 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
nd <- qcs_grid
nd$Y <- nd$Y * 10
expect_warning(p <- predict(fit, newdata = nd), regexp = "coordinates")
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.