Skip to content

Commit

Permalink
Fix #342
Browse files Browse the repository at this point in the history
Fix bug in `est_non_rf1` and `est_non_rf2`
columns when all the following conditions
were true:
  - predicting on new data
  - using a delta model
  - including IID random intercepts or
    time-varying coefficients
  • Loading branch information
seananderson committed May 17, 2024
1 parent 47ef042 commit a9bcda1
Show file tree
Hide file tree
Showing 4 changed files with 37 additions and 3 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Type: Package
Package: sdmTMB
Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
Version: 0.5.0.9006
Version: 0.5.0.9007
Authors@R: c(
person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca",
role = c("aut", "cre"),
Expand Down
7 changes: 7 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,12 @@
# sdmTMB (development version)

* Fix bug in `est_non_rf1` and `est_non_rf2` columns when all the following
conditions were true:
- predicting on new data
- using a delta model
- including IID random intercepts or time-varying coefficients
See #342. Thanks to @tom-peatman for the issue report.

* Add suggestion to use an optimized BLAS library to README.

* Add warning if it's detected that there were problems reloading (e.g., with
Expand Down
4 changes: 2 additions & 2 deletions src/sdmTMB.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1087,7 +1087,7 @@ Type objective_function<Type>::operator()()
if (!exclude_RE(k)) proj_iid_re_i(i,m) += RE(proj_RE_indexes(i, k) + temp,m);
}
}
proj_fe(i) += proj_iid_re_i(i,m);
proj_fe(i,m) += proj_iid_re_i(i,m);
}
}

Expand All @@ -1099,7 +1099,7 @@ Type objective_function<Type>::operator()()
for (int i = 0; i < proj_X_rw_ik.rows(); i++) {
for (int k = 0; k < proj_X_rw_ik.cols(); k++) {
proj_rw_i(i,m) += proj_X_rw_ik(i, k) * b_rw_t(proj_year(i), k, m);
proj_fe(i) += proj_rw_i(i,m);
proj_fe(i,m) += proj_rw_i(i,m);
}
}
}
Expand Down
27 changes: 27 additions & 0 deletions tests/testthat/test-random-intercepts.R
Original file line number Diff line number Diff line change
Expand Up @@ -297,3 +297,30 @@ test_that("Random intercept classes in predict() are checked appropriately", {
expect_s3_class(p11, "tbl_df")
})

test_that("Random intercepts are incorporated into est_non_rf1 and est_non_rf2 correctly #342", {
## test based on example by @tom-peatman #342
skip_on_cran()
skip_on_ci()
pcod_2011$year <- factor(pcod_2011$year)
## fit model with random intercepts for year in both model components
fit <- sdmTMB(
density ~ (1 | year),
spatial = "off",
data = pcod_2011, mesh = pcod_mesh_2011,
family = delta_gamma(link1 = "logit", link2 = "log")
)
## fitted random intercepts
t1 <- tidy(fit, effects = "ran_vals", model = 1)
t2 <- tidy(fit, effects = "ran_vals", model = 2)
expect_equal(t1$estimate, c(-0.24066, 0.33083, 0.20459, -0.29288), tolerance = 0.001)
expect_equal(t2$estimate, c(0.17724, -0.14142, 0.11862, -0.16844), tolerance = 0.001)

## data frame to get predictions for each year (at a specific location and depth)
new_dat <- expand.grid(
X = mean(pcod_2011$X), Y = mean(pcod_2011$Y),
depth = mean(pcod_2011$depth), year = unique(pcod_2011$year)
)
p <- predict(fit, newdata = new_dat)
expect_equal(p$est_non_rf1, c(-0.40011, 0.17137, 0.04514, -0.45234), tolerance = 0.001)
expect_equal(p$est_non_rf2, c(4.6455, 4.32684, 4.58688, 4.29982), tolerance = 0.001)
})

0 comments on commit a9bcda1

Please sign in to comment.