From 9d682b79ae318f475ee58860b4384f17cb604ba7 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 17 Oct 2023 22:15:46 -0700 Subject: [PATCH] Move add_barrier_mesh() to sdmTMBextra --- .github/workflows/R-CMD-check.yaml | 1 - .github/workflows/pkgdown.yaml | 2 +- DESCRIPTION | 5 - R/barrier.R | 162 ----------------------------- R/deprecated.R | 35 ++++++- R/fit.R | 12 +-- R/get-index-sims.R | 4 +- man/add_barrier_mesh.Rd | 103 ++---------------- man/dharma_residuals.Rd | 2 +- man/extract_mcmc.Rd | 2 +- man/get_index_sims.Rd | 4 +- man/sdmTMB.Rd | 12 +-- 12 files changed, 64 insertions(+), 280 deletions(-) delete mode 100644 R/barrier.R diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 32f578806..304910032 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -73,7 +73,6 @@ jobs: remotes::install_cran("sf") remotes::install_cran("rcmdcheck") remotes::install_cran("textshaping") - remotes::install_cran("rgdal") shell: Rscript {0} - name: Install dependencies Ubuntu diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index db8621735..3a2d9af98 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -31,7 +31,7 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::pkgdown, local::., any::rgdal, any::textshaping, any::remotes + extra-packages: any::pkgdown, local::., any::textshaping, any::remotes needs: website - name: Install extra packages for web-only vignettes diff --git a/DESCRIPTION b/DESCRIPTION index 3567a7a2f..b6ac52fdf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -92,13 +92,8 @@ Suggests: glmmTMB, ggplot2, inlabru, - INLAspacetime, - INLA, knitr, - rgdal, rmarkdown, - rnaturalearth, - rnaturalearthdata, sf, splancs, testthat, diff --git a/R/barrier.R b/R/barrier.R deleted file mode 100644 index 4654b20ba..000000000 --- a/R/barrier.R +++ /dev/null @@ -1,162 +0,0 @@ -#' Transform a mesh object into a mesh with correlation barriers -#' -#' @param spde_obj Output from [make_mesh()]. -#' @param barrier_sf An sf object with polygons defining the barriers. For -#' example, a coastline dataset for ocean data. **Note that this object must -#' have the same projection as the data used to generate the x and y columns -#' in `spde_obj`.** -#' @param range_fraction The fraction of the spatial range that barrier -#' triangles have. -#' @param proj_scaling If `spde_obj` was created with scaling of the coordinates -#' after the projection (e.g., dividing UTMs by 1000 so the spatial range is -#' on a reasonable scale) the x and y values in `spde_obj` are multiplied by -#' this scaling factor before applying the projection from `barrier_sf`. -#' @param plot Logical. -#' -#' @return A list similar to [make_mesh()] but with `spde_barrier` and a -#' couple other helper list elements added. -#' -#' If `plot = TRUE`, then a basic plot will be created as a side effect. Each -#' grey dot represents the center of a "normal" mesh triangle. Each red cross -#' represents the center of a "barrier" mesh triangle. -#' @export -#' @references -#' Bakka, H., Vanhatalo, J., Illian, J., Simpson, D., and Rue, H. 2019. -#' Non-stationary Gaussian models with physical barriers. -#' -#' -#' -#' -#' -#' @examples -#' if (require("sf", quietly = TRUE) && -#' require("ggplot2", quietly = TRUE) && -#' require("dplyr", quietly = TRUE) && -#' require("INLA", quietly = TRUE) && -#' require("rnaturalearth", quietly = TRUE)) { -#' -#' # For applied situations on finer scales, you may with to use scale = "large". -#' # For that, first: remotes::install_github("ropensci/rnaturalearthhires") -#' map_data <- rnaturalearth::ne_countries( -#' scale = "medium", -#' returnclass = "sf", country = "canada") -#' -#' # Crop the polygon for plotting and efficiency: -#' sf::st_bbox(map_data) -#' bc_coast <- suppressWarnings(suppressMessages( -#' sf::st_crop(map_data, -#' c(xmin = -134, ymin = 46, xmax = -120, ymax = 57)))) -#' -#' crs_utm9 <- 3156 # Pick a projection, here UTM9 -#' -#' sf::st_crs(bc_coast) <- 4326 # 'WGS84' -#' bc_coast <- sf::st_transform(bc_coast, crs_utm9) -#' -#' # Project our survey data coordinates: -#' survey <- pcod %>% dplyr::select(lon, lat, density) %>% -#' sf::st_as_sf(crs = 4326, coords = c("lon", "lat")) %>% -#' sf::st_transform(crs_utm9) -#' -#' # Plot our coast and survey data: -#' ggplot(bc_coast) + -#' geom_sf() + -#' geom_sf(data = survey, size = 0.5) -#' -#' # Note that a barrier mesh won't do much here for this -#' # example data set, but we nonetheless use it as an example. -#' -#' # Prepare for making the mesh -#' # First, we will extract the coordinates: -#' surv_utm_coords <- sf::st_coordinates(survey) -#' -#' # Then we will scale coordinates to km so the range parameter -#' # is on a reasonable scale for estimation: -#' pcod$X1000 <- surv_utm_coords[,1] / 1000 -#' pcod$Y1000 <- surv_utm_coords[,2] / 1000 -#' -# # Construct our mesh: -#' mesh <- make_mesh(pcod, xy_cols = c("X1000", "Y1000"), -#' n_knots = 200, type = "kmeans") -#' plot(mesh) -#' -#' # Add on the barrier mesh component: -#' bspde <- add_barrier_mesh( -#' mesh, bc_coast, range_fraction = 0.1, -#' proj_scaling = 1000, plot = TRUE -#' ) -#' -#' # In the above, the grey dots are the centre of triangles that are in the -#' # ocean. The red crosses are centres of triangles that are over land. The -#' # spatial range will be assumed to be 0.1 (`range_fraction`) over land -#' # compared to over water. -#' -#' # We can make a more advanced plot if we want: -#' mesh_df_water <- bspde$mesh_sf[bspde$normal_triangles, ] -#' mesh_df_land <- bspde$mesh_sf[bspde$barrier_triangles, ] -#' -#' ggplot(bc_coast) + -#' geom_sf() + -#' geom_sf(data = mesh_df_water, size = 1, colour = "blue") + -#' geom_sf(data = mesh_df_land, size = 1, colour = "green") -#' -#' # Now, when we fit our model with the new mesh, it will automatically -#' # include a barrier structure in the spatial correlation: -#' fit <- sdmTMB(density ~ s(depth, k = 3), data = pcod, mesh = bspde, -#' family = tweedie(link = "log")) -#' fit -#' } - -add_barrier_mesh <- function(spde_obj, barrier_sf, range_fraction = 0.2, - proj_scaling = 1, plot = FALSE) { - - if (!requireNamespace("sf", quietly = TRUE)) { - cli_abort("The sf package must be installed to use this function.") - } - if (!requireNamespace("INLAspacetime", quietly = TRUE)) { - cli_abort("The INLAspacetime package must be installed to use this function.") - } - if (!requireNamespace("INLAspacetime", quietly = TRUE)) { - cli_abort("The INLA package must be installed to use this function.") - } # b/c INLAspacetime uses it!! - - assert_that( - is.numeric(range_fraction), range_fraction <= 1, - range_fraction > 0, length(range_fraction) == 1 - ) - assert_that(is.numeric(proj_scaling), length(proj_scaling) == 1) - assert_that(inherits(barrier_sf, "sf")) - assert_that(inherits(spde_obj, "sdmTMBmesh")) - assert_that(is.logical(plot)) - - mesh <- spde_obj$mesh - tl <- length(mesh$graph$tv[, 1]) # the number of triangles in the mesh - pos_tri <- matrix(0, tl, 2) - for (i in seq_len(tl)) { - temp <- mesh$loc[mesh$graph$tv[i, ], ] - pos_tri[i, ] <- colMeans(temp)[c(1, 2)] - } - mesh_sf <- as.data.frame(pos_tri) - mesh_sf$X <- mesh_sf$V1 * proj_scaling - mesh_sf$Y <- mesh_sf$V2 * proj_scaling - mesh_sf <- sf::st_as_sf(mesh_sf, crs = sf::st_crs(barrier_sf), coords = c("X", "Y")) - - intersected <- sf::st_intersects(mesh_sf, barrier_sf) - water.triangles <- which(lengths(intersected) == 0) - land.triangles <- which(lengths(intersected) > 0) - - if (plot) { - plot(pos_tri[water.triangles, ], - col = "grey40", asp = 1, - xlab = spde_obj$xy_cols[1], ylab = spde_obj$xy_cols[2] - ) - points(pos_tri[land.triangles, ], col = "red", pch = 4) - } - # barrier_spde <- INLA::inla.barrier.fem(mesh, barrier.triangles = land.triangles) - barrier_spde <- INLAspacetime::mesh2fem.barrier(mesh, barrier.triangles = land.triangles) - spde_obj$spde_barrier <- barrier_spde - spde_obj$barrier_scaling <- c(1, range_fraction) - spde_obj$mesh_sf <- mesh_sf - spde_obj$barrier_triangles <- land.triangles - spde_obj$normal_triangles <- water.triangles - spde_obj -} diff --git a/R/deprecated.R b/R/deprecated.R index 8705bd7a0..456049e5c 100644 --- a/R/deprecated.R +++ b/R/deprecated.R @@ -3,7 +3,7 @@ #' Moved to the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} #' package #' -#' @param object **Deprecated** See the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package. +#' @param object **Deprecated** See the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package. Make sure to load \pkg{sdmTMBextra} *after* \pkg{sdmTMB}. #' #' @return #' Deprecated. See the @@ -22,7 +22,7 @@ extract_mcmc <- function(object = deprecated()) { #' DHARMa residuals #' #' Moved to the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} -#' package +#' package. Make sure to load \pkg{sdmTMBextra} *after* \pkg{sdmTMB}. #' #' @param simulated_response **Deprecated** See the \pkg{sdmTMBextra} package. #' @param object **Deprecated** See the \pkg{sdmTMBextra} package. @@ -42,3 +42,34 @@ extract_mcmc <- function(object = deprecated()) { dharma_residuals <- function(simulated_response = deprecated(), object, plot = TRUE, ...) { deprecate_stop("0.2.2", "dharma_residuals()", details = "Load the sdmTMBextra package to use this function: https://github.com/pbs-assess/sdmTMBextra") } + +#' Transform a mesh object into a mesh with correlation barriers +#' +#' Moved to the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} +#' package. Make sure to load \pkg{sdmTMBextra} *after* \pkg{sdmTMB}. +#' +#' @param spde_obj Output from [make_mesh()]. +#' @param barrier_sf An sf object with polygons defining the barriers. For +#' example, a coastline dataset for ocean data. **Note that this object must +#' have the same projection as the data used to generate the x and y columns +#' in `spde_obj`.** +#' @param range_fraction The fraction of the spatial range that barrier +#' triangles have. +#' @param proj_scaling If `spde_obj` was created with scaling of the coordinates +#' after the projection (e.g., dividing UTMs by 1000 so the spatial range is +#' on a reasonable scale) the x and y values in `spde_obj` are multiplied by +#' this scaling factor before applying the projection from `barrier_sf`. +#' @param plot Logical. +#' @export +#' @keywords internal +#' @return +#' Deprecated. See the +#' \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package. +#' @examples +#' \dontrun{ +#' add_barrier_mesh() +#' } +add_barrier_mesh <- function(spde_obj = deprecated(), barrier_sf = deprecated(), range_fraction = 0.2, + proj_scaling = 1, plot = FALSE) { + deprecate_stop("0.3.0", "add_barrier_mesh()", details = "Load the sdmTMBextra package to use this function: https://github.com/pbs-assess/sdmTMBextra \n We may move this back to the main package once the functionality joins fmesher. For now, this function requires INLA.") +} diff --git a/R/fit.R b/R/fit.R index 52e7ebfeb..85f222855 100644 --- a/R/fit.R +++ b/R/fit.R @@ -389,11 +389,6 @@ NULL #' # Perform several 'sanity' checks: #' sanity(fit) #' -#' # Visualize depth effect: (see ?visreg_delta) -#' visreg::visreg(fit, xvar = "depth") # link space; randomized quantile residuals -#' visreg::visreg(fit, xvar = "depth", scale = "response") -#' visreg::visreg(fit, xvar = "depth", scale = "response", gg = TRUE, rug = FALSE) -#' #' # Predict on the fitted data; see ?predict.sdmTMB #' p <- predict(fit) #' @@ -401,6 +396,12 @@ NULL #' p <- predict(fit, newdata = qcs_grid) #' head(p) #' +#' \donttest{ +#' # Visualize depth effect: (see ?visreg_delta) +#' visreg::visreg(fit, xvar = "depth") # link space; randomized quantile residuals +#' visreg::visreg(fit, xvar = "depth", scale = "response") +#' visreg::visreg(fit, xvar = "depth", scale = "response", gg = TRUE, rug = FALSE) +#' #' # Add spatiotemporal random fields: #' fit <- sdmTMB( #' density ~ 0 + as.factor(year), @@ -410,7 +411,6 @@ NULL #' ) #' fit #' -#' \donttest{ #' # Make the fields AR1: #' fit <- sdmTMB( #' density ~ s(depth), diff --git a/R/get-index-sims.R b/R/get-index-sims.R index 36be9657e..137ebb521 100644 --- a/R/get-index-sims.R +++ b/R/get-index-sims.R @@ -54,13 +54,14 @@ #' #' @export #' @examples -#' m <- sdmTMB(density ~ 0 + as.factor(year) + depth_scaled + depth_scaled2, +#' m <- sdmTMB(density ~ 0 + as.factor(year), #' data = pcod_2011, mesh = pcod_mesh_2011, family = tweedie(link = "log"), #' time = "year" #' ) #' qcs_grid_2011 <- replicate_df(qcs_grid, "year", unique(pcod_2011$year)) #' p <- predict(m, newdata = qcs_grid_2011, nsim = 100) #' x <- get_index_sims(p) +#' \donttest{ #' x_sims <- get_index_sims(p, return_sims = TRUE) #' #' if (require("ggplot2", quietly = TRUE)) { @@ -77,6 +78,7 @@ #' agg_function = function(x) sum(x), #' area_function = function(x, area) x * area #' ) +#' } get_index_sims <- function(obj, level = 0.95, return_sims = FALSE, diff --git a/man/add_barrier_mesh.Rd b/man/add_barrier_mesh.Rd index fc91ac56a..b69e86022 100644 --- a/man/add_barrier_mesh.Rd +++ b/man/add_barrier_mesh.Rd @@ -1,12 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/barrier.R +% Please edit documentation in R/deprecated.R \name{add_barrier_mesh} \alias{add_barrier_mesh} \title{Transform a mesh object into a mesh with correlation barriers} \usage{ add_barrier_mesh( - spde_obj, - barrier_sf, + spde_obj = deprecated(), + barrier_sf = deprecated(), range_fraction = 0.2, proj_scaling = 1, plot = FALSE @@ -31,99 +31,16 @@ this scaling factor before applying the projection from \code{barrier_sf}.} \item{plot}{Logical.} } \value{ -A list similar to \code{\link[=make_mesh]{make_mesh()}} but with \code{spde_barrier} and a -couple other helper list elements added. - -If \code{plot = TRUE}, then a basic plot will be created as a side effect. Each -grey dot represents the center of a "normal" mesh triangle. Each red cross -represents the center of a "barrier" mesh triangle. +Deprecated. See the +\href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package. } \description{ -Transform a mesh object into a mesh with correlation barriers +Moved to the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} +package. Make sure to load \pkg{sdmTMBextra} \emph{after} \pkg{sdmTMB}. } \examples{ -if (require("sf", quietly = TRUE) && - require("ggplot2", quietly = TRUE) && - require("dplyr", quietly = TRUE) && - require("INLA", quietly = TRUE) && - require("rnaturalearth", quietly = TRUE)) { - -# For applied situations on finer scales, you may with to use scale = "large". -# For that, first: remotes::install_github("ropensci/rnaturalearthhires") -map_data <- rnaturalearth::ne_countries( - scale = "medium", - returnclass = "sf", country = "canada") - -# Crop the polygon for plotting and efficiency: -sf::st_bbox(map_data) -bc_coast <- suppressWarnings(suppressMessages( - sf::st_crop(map_data, - c(xmin = -134, ymin = 46, xmax = -120, ymax = 57)))) - -crs_utm9 <- 3156 # Pick a projection, here UTM9 - -sf::st_crs(bc_coast) <- 4326 # 'WGS84' -bc_coast <- sf::st_transform(bc_coast, crs_utm9) - -# Project our survey data coordinates: -survey <- pcod \%>\% dplyr::select(lon, lat, density) \%>\% - sf::st_as_sf(crs = 4326, coords = c("lon", "lat")) \%>\% - sf::st_transform(crs_utm9) - -# Plot our coast and survey data: -ggplot(bc_coast) + - geom_sf() + - geom_sf(data = survey, size = 0.5) - -# Note that a barrier mesh won't do much here for this -# example data set, but we nonetheless use it as an example. - -# Prepare for making the mesh -# First, we will extract the coordinates: -surv_utm_coords <- sf::st_coordinates(survey) - -# Then we will scale coordinates to km so the range parameter -# is on a reasonable scale for estimation: -pcod$X1000 <- surv_utm_coords[,1] / 1000 -pcod$Y1000 <- surv_utm_coords[,2] / 1000 - -mesh <- make_mesh(pcod, xy_cols = c("X1000", "Y1000"), - n_knots = 200, type = "kmeans") -plot(mesh) - -# Add on the barrier mesh component: -bspde <- add_barrier_mesh( - mesh, bc_coast, range_fraction = 0.1, - proj_scaling = 1000, plot = TRUE -) - -# In the above, the grey dots are the centre of triangles that are in the -# ocean. The red crosses are centres of triangles that are over land. The -# spatial range will be assumed to be 0.1 (`range_fraction`) over land -# compared to over water. - -# We can make a more advanced plot if we want: -mesh_df_water <- bspde$mesh_sf[bspde$normal_triangles, ] -mesh_df_land <- bspde$mesh_sf[bspde$barrier_triangles, ] - -ggplot(bc_coast) + - geom_sf() + - geom_sf(data = mesh_df_water, size = 1, colour = "blue") + - geom_sf(data = mesh_df_land, size = 1, colour = "green") - -# Now, when we fit our model with the new mesh, it will automatically -# include a barrier structure in the spatial correlation: -fit <- sdmTMB(density ~ s(depth, k = 3), data = pcod, mesh = bspde, - family = tweedie(link = "log")) -fit -} +\dontrun{ +add_barrier_mesh() } -\references{ -Bakka, H., Vanhatalo, J., Illian, J., Simpson, D., and Rue, H. 2019. -Non-stationary Gaussian models with physical barriers. -\url{https://arxiv.org/abs/1608.03787} - -\url{https://sites.google.com/a/r-inla.org/www/barrier-model} - -\url{https://haakonbakkagit.github.io/btopic107.html} } +\keyword{internal} diff --git a/man/dharma_residuals.Rd b/man/dharma_residuals.Rd index 7af664865..88d489a87 100644 --- a/man/dharma_residuals.Rd +++ b/man/dharma_residuals.Rd @@ -21,7 +21,7 @@ Deprecated. See the } \description{ Moved to the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} -package +package. Make sure to load \pkg{sdmTMBextra} \emph{after} \pkg{sdmTMB}. } \examples{ \dontrun{ diff --git a/man/extract_mcmc.Rd b/man/extract_mcmc.Rd index c07e5ccd4..8fe4b9bfe 100644 --- a/man/extract_mcmc.Rd +++ b/man/extract_mcmc.Rd @@ -7,7 +7,7 @@ extract_mcmc(object = deprecated()) } \arguments{ -\item{object}{\strong{Deprecated} See the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package.} +\item{object}{\strong{Deprecated} See the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package. Make sure to load \pkg{sdmTMBextra} \emph{after} \pkg{sdmTMB}.} } \value{ Deprecated. See the diff --git a/man/get_index_sims.Rd b/man/get_index_sims.Rd index 1c27bed50..b7ddd0313 100644 --- a/man/get_index_sims.Rd +++ b/man/get_index_sims.Rd @@ -74,13 +74,14 @@ This function does nothing more than summarize and reshape the matrix of simulation draws into a data frame. } \examples{ -m <- sdmTMB(density ~ 0 + as.factor(year) + depth_scaled + depth_scaled2, +m <- sdmTMB(density ~ 0 + as.factor(year), data = pcod_2011, mesh = pcod_mesh_2011, family = tweedie(link = "log"), time = "year" ) qcs_grid_2011 <- replicate_df(qcs_grid, "year", unique(pcod_2011$year)) p <- predict(m, newdata = qcs_grid_2011, nsim = 100) x <- get_index_sims(p) +\donttest{ x_sims <- get_index_sims(p, return_sims = TRUE) if (require("ggplot2", quietly = TRUE)) { @@ -98,6 +99,7 @@ ind <- get_index_sims( area_function = function(x, area) x * area ) } +} \seealso{ \code{\link[=get_index]{get_index()}} } diff --git a/man/sdmTMB.Rd b/man/sdmTMB.Rd index 8109cfff9..1176c5c1e 100644 --- a/man/sdmTMB.Rd +++ b/man/sdmTMB.Rd @@ -349,11 +349,6 @@ tidy(fit, effects = "ran_par", conf.int = TRUE) # Perform several 'sanity' checks: sanity(fit) -# Visualize depth effect: (see ?visreg_delta) -visreg::visreg(fit, xvar = "depth") # link space; randomized quantile residuals -visreg::visreg(fit, xvar = "depth", scale = "response") -visreg::visreg(fit, xvar = "depth", scale = "response", gg = TRUE, rug = FALSE) - # Predict on the fitted data; see ?predict.sdmTMB p <- predict(fit) @@ -361,6 +356,12 @@ p <- predict(fit) p <- predict(fit, newdata = qcs_grid) head(p) +\donttest{ +# Visualize depth effect: (see ?visreg_delta) +visreg::visreg(fit, xvar = "depth") # link space; randomized quantile residuals +visreg::visreg(fit, xvar = "depth", scale = "response") +visreg::visreg(fit, xvar = "depth", scale = "response", gg = TRUE, rug = FALSE) + # Add spatiotemporal random fields: fit <- sdmTMB( density ~ 0 + as.factor(year), @@ -370,7 +371,6 @@ fit <- sdmTMB( ) fit -\donttest{ # Make the fields AR1: fit <- sdmTMB( density ~ s(depth),