From 3701635b77c0f64a6cf2c2818d6150252dc5e397 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Thu, 22 Feb 2024 12:25:55 -0800 Subject: [PATCH] Add extra_time tests and clean up fix --- DESCRIPTION | 4 +- R/fit.R | 6 +-- R/utils.R | 10 ++--- tests/testthat/test-extra-time.R | 66 ++++++++++++++++++++++++++++++++ 4 files changed, 75 insertions(+), 11 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index c7f9be64c..95b310ffc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: sdmTMB Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB' -Version: 0.4.2.9004 +Version: 0.4.2.9005 Authors@R: c( person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca", role = c("aut", "cre"), @@ -108,5 +108,5 @@ Config/testthat/parallel: true Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 SystemRequirements: GNU make diff --git a/R/fit.R b/R/fit.R index a2617152f..5fef87134 100644 --- a/R/fit.R +++ b/R/fit.R @@ -767,14 +767,14 @@ sdmTMB <- function( } if (!is.null(extra_time)) { # for forecasting or interpolating - data <- expand_time(df = data, time_slices = extra_time, time_column = time, weights = weights, offset = offset, upr = upr) + data <- expand_time(df = data, time_slices = extra_time, time_column = time, + weights = weights, offset = offset, upr = upr) offset <- data[["__sdmTMB_offset__"]] # expanded weights <- data[["__weight_sdmTMB__"]] # expanded upr <- data[["__dcens_upr__"]] # expanded - # data[["__dcens_upr__"]] <- NULL spde$loc_xy <- as.matrix(data[,spde$xy_cols,drop=FALSE]) spde$A_st <- fmesher::fm_basis(spde$mesh, loc = spde$loc_xy) - spde$sdm_spatial_id <- seq(1, nrow(data)) # FIXME + spde$sdm_spatial_id <- seq(1, nrow(data)) # FIXME? } check_irregalar_time(data, time, spatiotemporal, time_varying) diff --git a/R/utils.R b/R/utils.R index 5d1492af3..9a20a2699 100644 --- a/R/utils.R +++ b/R/utils.R @@ -229,13 +229,11 @@ parse_threshold_formula <- function(formula, thresh_type_short = "lin_thresh", } expand_time <- function(df, time_slices, time_column, weights, offset, upr) { - df[["__weight_sdmTMB__"]] <- ifelse(!is.null(weights), weights, 1) - df[["__sdmTMB_offset__"]] <- ifelse(!is.null(offset), offset, 0) - df[["__dcens_upr__"]] <- ifelse(!is.null(upr), upr, NA_real_) + df[["__weight_sdmTMB__"]] <- if (!is.null(weights)) weights else 1 + df[["__sdmTMB_offset__"]] <- if (!is.null(offset)) offset else 0 + df[["__dcens_upr__"]] <- if (!is.null(upr)) upr else NA_real_ fake_df <- df[1L, , drop = FALSE] - fake_df[["__weight_sdmTMB__"]] <- 0 - fake_df[["__sdmTMB_offset__"]] <- 0 - fake_df[["__dcens_upr__"]] <- NA_real_ + fake_df[["__weight_sdmTMB__"]] <- 0 # IMPORTANT: this turns off these data in the likelihood missing_years <- time_slices[!time_slices %in% df[[time_column]]] fake_df <- do.call("rbind", replicate(length(missing_years), fake_df, simplify = FALSE)) fake_df[[time_column]] <- missing_years diff --git a/tests/testthat/test-extra-time.R b/tests/testthat/test-extra-time.R index c3a584ede..37c8e51e8 100644 --- a/tests/testthat/test-extra-time.R +++ b/tests/testthat/test-extra-time.R @@ -110,3 +110,69 @@ test_that("extra_time, newdata, get_index() work", { expect_identical(ind6$year, c(2003, 2004, 2005, 2007, 2009, 2011, 2013, 2015, 2017)) expect_equal(ind3$est, ind6$est, tolerance = 0.1) }) +test_that("extra time does not affect estimation (without time series estimation)", { + # there was a bad bug at one point where the likelihood (via the weights) + # wasn't getting turned off for the extra time data! + skip_on_cran() + skip_on_ci() + # adding extra time at beginning or end + m <- sdmTMB(present ~ depth_scaled, + family = binomial(), + data = pcod_2011, spatial = "on", mesh = pcod_mesh_2011 + ) + m1 <- sdmTMB( + present ~ depth_scaled, + family = binomial(), data = pcod_2011, spatial = "on", + mesh = pcod_mesh_2011, + extra_time = 1990 + ) + m2 <- sdmTMB( + present ~ depth_scaled, + family = binomial(), data = pcod_2011, spatial = "on", + mesh = pcod_mesh_2011, + extra_time = 3000 + ) + expect_equal(m$model$par, m1$model$par) + expect_equal(m$model$par, m2$model$par) + + # with weights + set.seed(1) + w <- rlnorm(nrow(pcod_2011), meanlog = log(1), sdlog = 0.1) + m <- sdmTMB(present ~ depth_scaled, + family = binomial(), weights = w, + data = pcod_2011, spatial = "on", mesh = pcod_mesh_2011 + ) + m1 <- sdmTMB( + present ~ depth_scaled, weights = w, + family = binomial(), data = pcod_2011, spatial = "on", + mesh = pcod_mesh_2011, + extra_time = 1990 + ) + expect_equal(m$model$par, m1$model$par) + + # with offset as numeric + o <- log(w) + m <- sdmTMB(density ~ depth_scaled, + family = tweedie(), offset = o, + data = pcod_2011, mesh = pcod_mesh_2011 + ) + m1 <- sdmTMB(density ~ depth_scaled, + family = tweedie(), offset = o, + data = pcod_2011, mesh = pcod_mesh_2011, + extra_time = 1990 + ) + expect_equal(m$model$par, m1$model$par) + + # with offset as character + pcod_2011$off <- o + m <- sdmTMB(density ~ depth_scaled, + family = tweedie(), offset = "off", + data = pcod_2011, mesh = pcod_mesh_2011 + ) + m1 <- sdmTMB(density ~ depth_scaled, + family = tweedie(), offset = "off", + data = pcod_2011, mesh = pcod_mesh_2011, + extra_time = 1990 + ) + expect_equal(m$model$par, m1$model$par) +})