From 3167e0a14500b9e45bcd70ee1334f241e530bc05 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 3 Nov 2023 17:17:32 -0700 Subject: [PATCH] Fix offsets, extra time, and nsim combination #273 #274 --- R/predict.R | 11 +++- tests/testthat/test-extra-time.R | 12 ++++ tests/testthat/test-offset.R | 108 +++++++++++++++++++++++++------ 3 files changed, 110 insertions(+), 21 deletions(-) diff --git a/R/predict.R b/R/predict.R index c77d9613b..4714b40c3 100644 --- a/R/predict.R +++ b/R/predict.R @@ -312,8 +312,11 @@ predict.sdmTMB <- function(object, newdata = NULL, # places where we force newdata: nd_arg_was_null <- FALSE if (is.null(newdata)) { - if (is_delta(object) || nsim > 0 || type == "response" || !is.null(mcmc_samples) || se_fit || !is.null(re_form) || !is.null(re_form_iid)) { + if (is_delta(object) || nsim > 0 || type == "response" || !is.null(mcmc_samples) || se_fit || !is.null(re_form) || !is.null(re_form_iid) || !is.null(offset)) { newdata <- object$data + if (!is.null(object$extra_time)) { # issue #273 + newdata <- newdata[!newdata[[object$time]] %in% object$extra_time,] + } nd_arg_was_null <- TRUE # will be used to carry over the offset } } @@ -341,6 +344,7 @@ predict.sdmTMB <- function(object, newdata = NULL, tmb_data <- object$tmb_data tmb_data$do_predict <- 1L no_spatial <- as.logical(object$tmb_data$no_spatial) + fake_nd <- NULL if (!is.null(newdata)) { if (any(!xy_cols %in% names(newdata)) && isFALSE(pop_pred) && !no_spatial) @@ -494,7 +498,7 @@ predict.sdmTMB <- function(object, newdata = NULL, cli_abort("Prediction offset vector does not equal number of rows in prediction dataset.") } tmb_data$proj_offset_i <- if (!is.null(offset)) offset else rep(0, nrow(proj_X_ij[[1]])) - if (nd_arg_was_null) tmb_data$proj_offset_i <- tmb_data$offset_i + # if (nd_arg_was_null) tmb_data$proj_offset_i <- tmb_data$offset_i tmb_data$proj_X_threshold <- thresh[[1]]$X_threshold # TODO DELTA HARDCODED TO 1 tmb_data$area_i <- if (length(area) == 1L) rep(area, nrow(proj_X_ij[[1]])) else area tmb_data$proj_mesh <- proj_mesh @@ -664,6 +668,9 @@ predict.sdmTMB <- function(object, newdata = NULL, } } + if (!is.null(fake_nd)) { + out <- out[-seq(nrow(out) - nrow(fake_nd) + 1, nrow(out)), ,drop=FALSE] # issue #273 + } return(out) } diff --git a/tests/testthat/test-extra-time.R b/tests/testthat/test-extra-time.R index 0c5f27062..431ec70ee 100644 --- a/tests/testthat/test-extra-time.R +++ b/tests/testthat/test-extra-time.R @@ -24,4 +24,16 @@ test_that("extra time, newdata, and offsets work", { expect_equal(nrow(p4), nrow(pcod)) expect_equal(p1$est, p2$est) expect_equal(p3$est, p4$est) + + #273 (with nsim) + set.seed(1) + suppressWarnings(p5 <- predict(m, offset = pcod$os, nsim = 2L)) + expect_equal(ncol(p5), 2L) + expect_equal(nrow(p5), nrow(pcod)) + + set.seed(1) + suppressWarnings(p6 <- predict(m, newdata = pcod, offset = pcod$os, nsim = 2L)) + expect_equal(ncol(p6), 2L) + expect_equal(nrow(p6), nrow(pcod)) + expect_equal(p6[,1,drop=TRUE], p5[,1,drop=TRUE]) }) diff --git a/tests/testthat/test-offset.R b/tests/testthat/test-offset.R index 2a5ad8332..ddeaf5fae 100644 --- a/tests/testthat/test-offset.R +++ b/tests/testthat/test-offset.R @@ -29,15 +29,13 @@ test_that("Offset matches glmmTMB", { family = Gamma(link = "log") ) - fit1_off <- glmmTMB::glmmTMB(density ~ 1, - offset = pcod_pos$offset, + fit1_off <- glmmTMB::glmmTMB(density ~ 1 + offset(offset), data = pcod_pos, family = Gamma(link = "log") ) pcod_pos$offset2 <- log(1) - fit1_off0 <- glmmTMB::glmmTMB(density ~ 1, - offset = pcod_pos$offset2, + fit1_off0 <- glmmTMB::glmmTMB(density ~ 1 + offset(offset2), data = pcod_pos, family = Gamma(link = "log") ) @@ -103,41 +101,113 @@ test_that("Offset prediction matches glm()", { skip_on_cran() skip_if_not_installed("glmmTMB") set.seed(1) - pcod$offset <- rnorm(nrow(pcod)) + dat <- pcod[pcod$density > 0,] + dat$.offset <- rnorm(nrow(dat)) fit <- sdmTMB( present ~ 1, - offset = pcod$offset, - data = pcod, spatial = "off", - family = binomial() + offset = dat$.offset, + data = dat, spatial = "off", + family = Gamma(link = "log") ) fit_glm <- glm( - present ~ 1, - offset = pcod$offset, - data = pcod, - family = binomial() + present ~ 1 + offset(.offset), + data = dat, + family = Gamma(link = "log") ) fit_glmmTMB <- glmmTMB::glmmTMB( - present ~ 1, - offset = pcod$offset, - data = pcod, - family = binomial() + present ~ 1 + offset(.offset), + data = dat, + family = Gamma(link = "log") ) p <- predict(fit) p_glm <- predict(fit_glm) p_glmmTMB <- predict(fit_glmmTMB) + expect_equal(p$est, unname(p_glm)) + expect_equal(p$est, p_glmmTMB) + p_glmmTMB <- predict(fit_glmmTMB, newdata = dat) expect_equal(p$est, unname(p_glm)) expect_equal(p$est, p_glmmTMB) set.seed(1) - p <- predict(fit, nsim = 1000) + p <- predict(fit, nsim = 1000, offset = dat$.offset) mu <- apply(p, 1, mean) plot(mu, p_glm) expect_equal(unname(mu), unname(p_glm), tolerance = 0.01) # sdmTMB ignores offset here (but not glm() or glmmTMB()!) - # p <- predict(fit, newdata = pcod) - # p_glmmTMB <- predict(fit_glmmTMB, newdata = pcod) + # p <- predict(fit, newdata = dat) + # p_glmmTMB <- predict(fit_glmmTMB, newdata = dat) # expect_equal(p$est, unname(p_glmmTMB)) }) +# # +# # offset/prediction setting checks: +# +# pos <- dogfish[dogfish$catch_weight > 0,] +# +# m1 <- sdmTMB(catch_weight ~ 1, family = Gamma("log"), data = pos, offset = log(pos$area_swept), spatial = "off") +# # m2 <- glmmTMB::glmmTMB(catch_weight ~ 1, family = Gamma("log"), data = pos, offset = log(pos$area_swept)) +# # m3 <- glm(catch_weight ~ 1, family = Gamma("log"), data = pos, offset = log(pos$area_swept)) +# +# head(predict(m1, newdata = pos, offset = rep(0, nrow(pos)))$est) # right +# head(predict(m1, offset = rep(0, nrow(pos)))$est) # right (was wrong) +# +# head(predict(m1, offset = log(pos$area_swept))$est) # right +# head(predict(m1, newdata = pos, offset = log(pos$area_swept))$est) # right +# +# head(predict(m1)$est) # right +# head(predict(m1, newdata = pos)$est) # right +# +# head(predict(m1, newdata = pos, offset = rep(0, nrow(pos)), nsim = 2)) # right +# head(predict(m1, offset = rep(0, nrow(pos)), nsim = 2)) # right (was wrong) +# +# head(predict(m1, offset = log(pos$area_swept), nsim = 2)) # right +# head(predict(m1, newdata = pos, offset = log(pos$area_swept), nsim = 2)) # right +# +# head(predict(m1, nsim = 2)) # right +# head(predict(m1, newdata = pos, nsim = 2)) # right +# +# +# # m2 <- glmmTMB::glmmTMB(catch_weight ~ 1, family = Gamma("log"), data = pos) +# # head(predict(m2)) +# # head(predict(m2, newdata = pos)) +# # predict(m2, newdata = pos[1:3,,drop=FALSE]) +# # +# # e1 <- predict(m1, newdata = pos, offset = rep(0, nrow(pos)))$est +# # plot(exp(e1), pos$catch_weight) +# # mean(exp(e1)) +# # mean(pos$catch_weight) +# # +# # +# # e1 <- predict(m1, newdata = pos, offset = log(pos$area_swept))$est +# # plot(exp(e1), pos$catch_weight) +# # mean(exp(e1)) +# # mean(pos$catch_weight) +# # +# # head(predict(m1, newdata = pos, offset = rep(0, nrow(pos)))$est) +# # +# # head(predict(m2)) +# # head(predict(m2, newdata = pos)) +# # +# # head(predict(m2)) +# # head(predict(m2, newdata = pos)) +# # head(predict(m3)) +# # +# # expect_equal(p7[,1,drop=TRUE], p6[,1,drop=TRUE]) +# # +# # set.seed(1) +# # suppressWarnings(p8 <- predict(m, newdata = pcod, nsim = 2L)) +# # suppressWarnings(p9 <- predict(m, newdata = pcod, offset = rep(0, nrow(pcod)), nsim = 2L)) +# # +# # pos <- subset(dogfish, catch_weight > 0) +# # +# # mm <- glm(catch_weight ~ 1, family = Gamma("log"), data = pos, offset = log(pos$area_swept)) +# # mm_s <- sdmTMB(catch_weight ~ 1, family = Gamma("log"), data = pos, offset = log(pos$area_swept), spatial = "off") +# # pp1 <- predict(mm) +# # pp1_s <- predict(mm_s) +# # +# # +# # +# # head(p8) +# # head(p9)