diff --git a/DESCRIPTION b/DESCRIPTION index 5a78d950..30271b6f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: sdmTMB Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB' -Version: 0.6.0.9013 +Version: 0.6.0.9015 Authors@R: c( person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca", role = c("aut", "cre"), @@ -64,7 +64,6 @@ Depends: Imports: assertthat, abind, - clisymbols, cli, fmesher, fishMod, diff --git a/NEWS.md b/NEWS.md index 27f53b2a..ef4c599a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,13 @@ # sdmTMB (development version) +* Fix bug in `est` column when predicting on new data with Poisson-link + delta models with `type = "link"` and `re_form = NA`. #389 + +* Fix bug in `s95` param reporting from the `tidy()` method. `s95` is present + in the logistic threshold models. The model itself was fine but the `s95` + parameter was supposed to be reported by `tidy()` as a combination of two + other parameters. This also affected the output in `print()`/`summary()`. + * Add `return_tmb_report` to `simulate.sdmTMB()`. * Add `newdata` argument to `simulate.sdmTMB()`. This enables simulating on diff --git a/R/mesh.R b/R/mesh.R index 0a071ed1..4b6086eb 100644 --- a/R/mesh.R +++ b/R/mesh.R @@ -207,12 +207,12 @@ binary_search_knots <- function(loc_xy, cat("cutoff =", pretty_cutoff, "| knots =", realized_knots) if (realized_knots > n_knots) { L <- m + 1 - cat(" |", clisymbols::symbol$arrow_down, "\n") + cat(" |", cli::symbol$arrow_down, "\n") } else if (realized_knots < n_knots) { R <- m - 1 - cat(" |", clisymbols::symbol$arrow_up, "\n") + cat(" |", cli::symbol$arrow_up, "\n") } else { - cat(" |", clisymbols::symbol$tick, "\n") + cat(" |", cli::symbol$tick, "\n") return(mesh) } } diff --git a/R/predict.R b/R/predict.R index dac6a165..d3b0aab3 100644 --- a/R/predict.R +++ b/R/predict.R @@ -785,14 +785,28 @@ predict.sdmTMB <- function(object, newdata = NULL, if (type == "response") { nd$est1 <- object$family[[1]]$linkinv(r$proj_fe[,1]) nd$est2 <- object$family[[2]]$linkinv(r$proj_fe[,2]) - nd$est <- nd$est1 * nd$est2 + if (object$tmb_data$poisson_link_delta) { + .n <- nd$est1 # expected group density (already exp()) + .p <- 1 - exp(-.n) # expected encounter rate + .w <- nd$est2 # expected biomass per group (already exp()) + .r <- (.n * .w) / .p # (n * w)/p # positive expectation + nd$est1 <- .p # expected encounter rate + nd$est2 <- .r # positive expectation + nd$est <- .n * .w # expected combined value + } else { + nd$est <- nd$est1 * nd$est2 + } } else { nd$est1 <- r$proj_fe[,1] nd$est2 <- r$proj_fe[,2] if (is.na(model)) { p1 <- object$family[[1]]$linkinv(r$proj_fe[,1]) p2 <- object$family[[2]]$linkinv(r$proj_fe[,2]) - nd$est <- object$family[[2]]$linkfun(p1 * p1) + if (object$tmb_data$poisson_link_delta) { + nd$est <- nd$est1 + nd$est2 + } else { + nd$est <- object$family[[2]]$linkfun(p1 * p1) + } if (se_fit) { nd$est <- sr_est_rep$proj_rf_delta nd$est_se <- sr_se_rep$proj_rf_delta diff --git a/R/tidy.R b/R/tidy.R index 71fe69e0..257905ca 100644 --- a/R/tidy.R +++ b/R/tidy.R @@ -135,14 +135,18 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model = if (x$tmb_data$threshold_func > 0) { if (x$threshold_function == 1L) { par_name <- paste0(x$threshold_parameter, c("-slope", "-breakpt")) + estimates <- est$b_threshold[,model,drop=TRUE] + ses <- se$b_threshold[,model,drop=TRUE] } else { par_name <- paste0(x$threshold_parameter, c("-s50", "-s95", "-smax")) + estimates <- c(est$s50[model], est$s95[model], est$s_max[model]) + ses <- c(se$s50[model], se$s95[model], se$s_max[model]) } out <- rbind( out, data.frame( - term = par_name, estimate = est$b_threshold[,model,drop=TRUE], - std.error = se$b_threshold[,model,drop=TRUE], stringsAsFactors = FALSE + term = par_name, estimate = estimates, + std.error = ses, stringsAsFactors = FALSE ) ) } diff --git a/README.md b/README.md index 7147a42e..8e893844 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ -# sdmTMB +# sdmTMB > Spatial and spatiotemporal GLMMs with TMB diff --git a/header.md b/header.md index a2499511..d6766460 100644 --- a/header.md +++ b/header.md @@ -1,6 +1,6 @@ -# sdmTMB +# sdmTMB > Spatial and spatiotemporal GLMMs with TMB diff --git a/man/figures/logo-sdmTMB.png b/man/figures/logo.png similarity index 100% rename from man/figures/logo-sdmTMB.png rename to man/figures/logo.png diff --git a/pkgdown/favicon/apple-touch-icon.png b/pkgdown/favicon/apple-touch-icon.png new file mode 100644 index 00000000..6722bda1 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon.png differ diff --git a/pkgdown/favicon/favicon-96x96.png b/pkgdown/favicon/favicon-96x96.png new file mode 100644 index 00000000..dfc95836 Binary files /dev/null and b/pkgdown/favicon/favicon-96x96.png differ diff --git a/pkgdown/favicon/favicon.ico b/pkgdown/favicon/favicon.ico new file mode 100644 index 00000000..f0280701 Binary files /dev/null and b/pkgdown/favicon/favicon.ico differ diff --git a/pkgdown/favicon/favicon.svg b/pkgdown/favicon/favicon.svg new file mode 100644 index 00000000..23e20208 --- /dev/null +++ b/pkgdown/favicon/favicon.svg @@ -0,0 +1,3 @@ + \ No newline at end of file diff --git a/pkgdown/favicon/site.webmanifest b/pkgdown/favicon/site.webmanifest new file mode 100644 index 00000000..4ebda26b --- /dev/null +++ b/pkgdown/favicon/site.webmanifest @@ -0,0 +1,21 @@ +{ + "name": "", + "short_name": "", + "icons": [ + { + "src": "/web-app-manifest-192x192.png", + "sizes": "192x192", + "type": "image/png", + "purpose": "maskable" + }, + { + "src": "/web-app-manifest-512x512.png", + "sizes": "512x512", + "type": "image/png", + "purpose": "maskable" + } + ], + "theme_color": "#ffffff", + "background_color": "#ffffff", + "display": "standalone" +} \ No newline at end of file diff --git a/pkgdown/favicon/web-app-manifest-192x192.png b/pkgdown/favicon/web-app-manifest-192x192.png new file mode 100644 index 00000000..47b1e6b5 Binary files /dev/null and b/pkgdown/favicon/web-app-manifest-192x192.png differ diff --git a/pkgdown/favicon/web-app-manifest-512x512.png b/pkgdown/favicon/web-app-manifest-512x512.png new file mode 100644 index 00000000..18f90b4f Binary files /dev/null and b/pkgdown/favicon/web-app-manifest-512x512.png differ diff --git a/scratch/binary-search.R b/scratch/binary-search.R index f53ed20d..64ebf678 100644 --- a/scratch/binary-search.R +++ b/scratch/binary-search.R @@ -47,12 +47,12 @@ binary_search_knots <- function(x, y, cat("cutoff =", pretty_cutoff, "| knots =", realized_knots) if (realized_knots > n_knots) { L <- m + 1 - cat(" |", clisymbols::symbol$arrow_down, "\n") + cat(" |", cli::symbol$arrow_down, "\n") } else if (realized_knots < n_knots) { R <- m - 1 - cat(" |", clisymbols::symbol$arrow_up, "\n") + cat(" |", cli::symbol$arrow_up, "\n") } else { - cat(" |", clisymbols::symbol$tick, "\n") + cat(" |", cli::symbol$tick, "\n") return(mesh) } } diff --git a/tests/testthat/test-poisson-link.R b/tests/testthat/test-poisson-link.R new file mode 100644 index 00000000..932047c9 --- /dev/null +++ b/tests/testthat/test-poisson-link.R @@ -0,0 +1,23 @@ +test_that("Poisson-link prediction issues from GitHub issues #389,", { + skip_on_cran() + dogfish$log_depth <- log(dogfish$depth) + mesh <- make_mesh(dogfish, c("X", "Y"), cutoff = 30) + fit_dpg <- sdmTMB(catch_weight ~ 0 + as.factor(year) + s(log_depth), + family = delta_gamma(type = "poisson-link"), + spatial = "on", + mesh = mesh, + data = dogfish, + offset = log(dogfish$area_swept) + ) + p <- predict(fit_dpg, re_form = NA) + p$est_est1est2 <- p$est1 + p$est2 + expect_equal(p$est, p$est_est1est2) + expect_equal(p$est_est1est2[1:5], + c(7.01028301430728, 2.34143881755487, 6.96979232578834, 6.99973559970208, + 7.03187981132451), tolerance = 1e-3) + + pp <- predict(fit_dpg, type = "response") + expect_equal(pp$est[1:5], + c(693.502617933497, 107.803298324919, 8414.54507288536, 5770.52404422525, + 6545.06096568627), tolerance = 1e-3) +}) diff --git a/tests/testthat/test-threshold-models.R b/tests/testthat/test-threshold-models.R index db79cef5..15dfcb5d 100644 --- a/tests/testthat/test-threshold-models.R +++ b/tests/testthat/test-threshold-models.R @@ -10,7 +10,7 @@ test_that("A logistic threshold model fits", { expect_true("depth_scaled-s50" %in% tidy(m)$term) expect_true("depth_scaled-s95" %in% tidy(m)$term) expect_true("depth_scaled-smax" %in% tidy(m)$term) - expect_equal(tidy(m)[,"estimate",drop=TRUE], c(1.555 , 1.655 , 1.718 , 1.138, -0.979, -3.173 , 1.760), tolerance = 1e-3) + expect_equal(tidy(m)[,"estimate",drop=TRUE], c(1.555 , 1.655 , 1.718 , 1.138, -0.979, -0.937 , 1.760), tolerance = 1e-3) }) test_that("A linear threshold model fits", { @@ -114,3 +114,4 @@ test_that("A linear threshold *delta* model fits", { expect_equal(t1$std.error, td1$std.error, tolerance = 1e-5) expect_equal(t2$std.error, td2$std.error, tolerance = 1e-5) }) + diff --git a/vignettes/articles/ggeffects.Rmd b/vignettes/articles/ggeffects.Rmd index 72ee35ee..dc634117 100644 --- a/vignettes/articles/ggeffects.Rmd +++ b/vignettes/articles/ggeffects.Rmd @@ -81,7 +81,7 @@ plot(g2) We can add in data points ```{r} -plot(g, add.data = TRUE) +plot(g, show_data = TRUE) ``` We can also use `ggeffect` to plot multiple variables by listing them in `terms = c()`, with the first term listed indicating the variable to be plotted on the x-axis, and the remaining listed terms (up to four total) indicating the groups. @@ -120,7 +120,7 @@ g2 <- ggeffect(fit2, "depth_scaled [-3:2.7 by=0.05]") plot(g2) # note the high density values dwarf the fitted curve here -plot(g2, add.dat = TRUE) +plot(g2, show_data = TRUE) ``` We can fit a model with an interaction of two continuous variables: diff --git a/vignettes/web_only/multispecies.Rmd b/vignettes/web_only/multispecies.Rmd new file mode 100644 index 00000000..c61f4e33 --- /dev/null +++ b/vignettes/web_only/multispecies.Rmd @@ -0,0 +1,245 @@ +--- +title: "Fitting multispecies models with sdmTMB" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Fitting multispecies models with sdmTMB} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +**If the code in this vignette has not been evaluated, a rendered version is available on the [documentation site](https://pbs-assess.github.io/sdmTMB/index.html) under 'Articles'.** + +```{r setup, include = FALSE, cache = FALSE} +EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") +knitr::opts_chunk$set( + collapse = TRUE, + warning = FALSE, + message = FALSE, + comment = "#>", + fig.width = 7, + fig.asp = 0.618, + eval = EVAL, + purl = EVAL +) +``` + +```{r packages, message=FALSE} +library(sdmTMB) +``` + +For some applications, we might be interested in fitting a model that includes multiple responses such as 2+ species, or multiple size or age classes within a species. The most important step in fitting these models is understanding which parameters are shared, and which parameters are species-specific. + +Below, we illustrate a series of models. We'll start by simulating a 2-species dataset. Each species is allowed to have unique spatial standard deviations (`sigma_O`) as well as different year effects. + +```{r sim_dat} +set.seed(1) +predictor_dat <- data.frame( + X = runif(1000), Y = runif(1000), + year = rep(1:5, each = 200) +) +predictor_dat$fyear <- as.factor(predictor_dat$year) +mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1) +sim_dat_A <- sdmTMB_simulate( + formula = ~ 0 + fyear, + data = predictor_dat, + time = "year", + mesh = mesh, + range = 0.2, + family = gaussian(), + seed = 42, + sigma_O = 0.2, + phi = 0.1, + sigma_E = 0.3, + B = runif(5, min = 5, max = 8) # 5 random year effects +) +sim_dat_A$species <- "A" +sim_dat_B <- sdmTMB_simulate( + formula = ~ 0 + fyear, + data = predictor_dat, + time = "year", + mesh = mesh, + range = 0.2, + family = gaussian(), + seed = 43, + sigma_O = 0.3, + phi = 0.1, + sigma_E = 0.3, + B = runif(5, min = 5, max = 8) # 5 random year effects +) +sim_dat_B$species <- "B" +sim_dat <- rbind(sim_dat_A, sim_dat_B) +sim_dat$fyear <- factor(sim_dat$year) +``` + +We'll start by making an SPDE mesh across the full dataset. + +```{r mesh_fig, fig.asp=1} +mesh <- make_mesh(sim_dat, c("X", "Y"), cutoff = 0.1) +plot(mesh) +``` + +### Model 1: species-specific intercepts + +As a first model, we can include species-specific year effects. This can be done in a couple ways. One option would be to estimate the `species * year` interaction, letting the year effects for each species be independent. Here, all other parameters and random effect values (range, spatial field, spatial variance, spatiotemporal fields, spatiotemporal variances) are shared. + +```{r} +fit <- sdmTMB( + observed ~ fyear * species, + data = sim_dat, + time = "year", + spatiotemporal = "iid", + mesh = mesh, + family = gaussian() +) +fit +``` + +### Model 2: species-specific spatial fields + +We may be interested in fitting a model that lets the spatial patterning differ by species. These kinds of models can be expressed using spatially varying coefficients. Note that we use `spatial = off` because this represents a global spatial intercept---turning this off is akin to using a `-1` of `0` in a main formula including a factor. Both species take their spatial fields from the `spatial_varying` field here. + +```{r} +fit <- sdmTMB( + observed ~ fyear * species, + data = sim_dat, + mesh = mesh, + family = gaussian(), + spatial = "off", + time = "year", + spatiotemporal = "iid", + spatial_varying = ~ 0 + factor(species) +) +fit +``` + +You'll notice that there are two rows of entries for `sigma_Z` our spatially varying random field standard deviation: + +```{r} +tidy(fit, "ran_pars") +``` + +This means that our model is trying to estimate separate species-specific variance terms for the species-specific spatial fields (say *that* 10 times fast!). Here, that matches how we simulated the data. In other contexts, especially if we ran into estimation issues, we might want to share those SDs. + +If we wanted to estimate species-specific spatial fields with a single shared variance (meaning the net magnitude of the peaks and valleys in the fields were similar but the wiggles themselves were species specific), we could do that by specifying a custom map argument and passing it into `sdmTMBcontrol()`. Any shared factor levels in the `map` are gathered to have 'mirrored' or shared parameter values. We assign these to `ln_tau_Z` because, internally, this is the parameter that gets converted into the spatially-varying field variances (the SDs of those fields are `sigma_Z`). + +This case is pretty simple, but for more complicated cases we could figure out the structure of our needed `map` vector as follows: + +```{r} +colnames(model.matrix(~ 0 + factor(species), data = sim_dat)) +``` + +So, we need a vector of length two with shared factor values: + +```{r} +map_list <- list(ln_tau_Z = factor(c(1, 1))) +``` + +Then, we can use our map list to share the spatially varying coefficient SDs: + +```{r} +fit <- sdmTMB( + observed ~ fyear * factor(species), + data = sim_dat, + mesh = mesh, + family = gaussian(), + spatial = "off", + time = "year", + spatiotemporal = "iid", + spatial_varying = ~ 0 + factor(species), + control = sdmTMBcontrol(map = map_list) +) +fit +``` + +Notice the spatially varying coefficient SD is now shared. + +### Model 3: species-specific spatiotemporal fields + +In all of the examples above, spatiotemporal fields are included, but shared among species. As another example, we can extend the above approaches to set up a model that includes spatiotemporal fields unique to each species. + +One approach to including separate spatiotemporal fields by species is creating a new variable that is a concatenation of species and year (or any given time step factor). For example, we can then implement this form of species-specific spatiotemporal variation by changing the `time` argument to be `time = "species_year"`. + +```{r} +sim_dat$species_year <- factor(paste(sim_dat$species, sim_dat$year)) +map_list <- list(ln_tau_Z = factor(c(1, 1))) +fit <- sdmTMB( + observed ~ fyear * factor(species), + data = sim_dat, + mesh = mesh, + family = gaussian(), + spatial = "off", + time = "species_year", + spatiotemporal = "iid", + spatial_varying = ~ 0 + factor(species), + control = sdmTMBcontrol(map = map_list) +) +fit +``` + +### Model 4: hack species into the time element for spatial models + +If we only wanted to fit a spatial model but had several species (or other groups), one approach is to pretend our species (or other group) is the time element. + +```{r} +sim_dat$numeric_species <- as.numeric(factor(sim_dat$species)) # needs to be numeric +fit_fake_time <- sdmTMB( + observed ~ 0 + factor(species), + data = sim_dat, + mesh = mesh, + family = gaussian(), + spatial = "off", + time = "numeric_species", #< hack + spatiotemporal = "iid" #< 'AR1' or 'RW' probably wouldn't make sense here +) +fit_fake_time +``` + +This is just a convenience though. We could instead do the same thing using the `spatial_varying` argument making sure to 'map' the field variances to be shared to match the above: + +```{r} +fit_svc <- sdmTMB( + observed ~ 0 + factor(species), + data = sim_dat, + mesh = mesh, + family = gaussian(), + spatial = "off", + spatial_varying = ~ 0 + factor(species), + control = sdmTMBcontrol(map = list(ln_tau_Z = factor(c(1, 1)))) +) +fit_svc +``` + +We can prove they're identical: + +```{r} +logLik(fit_fake_time) +logLik(fit_svc) +``` + +### Putting it all together + +These examples illustrate a number of ways that species-specific effects can be included in `sdmTMB` models, and can be extended to other categories/groups/cohorts within a species for which one wants to control the amount of information shared between groups (e.g., age-, size-, or stage-specific estimates). A brief summary of these approaches can be summarized as: + +```{r echo=FALSE} +desc <- data.frame( + Form = c( + "Main effects", + "Spatial effects", + "Spatial effects w/shared variance", + "Spatiotemporal effects"), + Implementation = c( + "Year-by-species interactions or smooths by year", + "Spatially varying coefficients", + "Spatially varying coefficients + map argument", + "Species-year factor as time variable") +) +if (require("knitr", quietly = TRUE)) { + knitr::kable(desc) +} else + print(desc) +``` + +### Further extensions + +As long as you're willing to treat spatiotemporal and group-level fields (e.g., for different species or age cohorts) as independent, sdmTMB can be used to fit models to these data. For example, this allows sdmTMB to be used for standardization of age or length composition data as in [Thorson and Haltuch (2018) CJFAS](https://doi.org/10.1139/cjfas-2018-0015). The approach is to similar to the above and we plan to write a separate vignette on the topic.