diff --git a/DESCRIPTION b/DESCRIPTION index e3fad074..eaf29de5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -66,7 +66,8 @@ Suggests: glmnetUtils, gnlm, geosphere, - inlabru (>= 2.6.0), + inlabru (>= 2.10.0), + fmesher (>= 0.1.7), igraph, knitr, mboost, diff --git a/NEWS.md b/NEWS.md index e0b67e64..d7fff781 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ * :fire: Enable stars and multi-temporal SpatRaster zones for `scenario()` and `distribution()` #121 #### Minor improvements and bug fixes -* :bug: fix for support of factor x continious variable interaction #131 +* :bug: fix for support of factor x continuous variable interaction #131 * Renamed `add_control_extrapolation` to `add_limits_extrapolation()`. * :bug: fix to `engine_gdb` also to support non-linear smooth functions (again). * Small fix to support deprecated `field_occurrence` field in `validate` for convenience. diff --git a/R/add_latent.R b/R/add_latent.R index f2736f49..08b90241 100644 --- a/R/add_latent.R +++ b/R/add_latent.R @@ -84,7 +84,7 @@ methods::setMethod( methods::signature(x = "BiodiversityDistribution"), function(x, method = 'spde', priors = NULL, separate_spde = FALSE, ...) { assertthat::assert_that(inherits(x, "BiodiversityDistribution"), - is.character(method) && !missing(method), + is.character(method), is.null(priors) || inherits(priors, 'PriorList'), is.logical(separate_spde) ) diff --git a/R/add_predictors.R b/R/add_predictors.R index 0af345a0..2c227a13 100644 --- a/R/add_predictors.R +++ b/R/add_predictors.R @@ -164,7 +164,7 @@ methods::setMethod( } # Check that all names allowed - problematic_names <- grep("offset|w|weight|spatial_offset|Intercept|spatial.field", names(env),fixed = TRUE) + problematic_names <- grep("offset|w|weight|spatial_offset|observed|Intercept|spatial.field", names(env),fixed = TRUE) if( length(problematic_names)>0 ){ stop(paste0("Some predictor names are not allowed as they might interfere with model fitting:", paste0(names(env)[problematic_names],collapse = " | "))) } diff --git a/R/class-distributionmodel.R b/R/class-distributionmodel.R index 04f36aa7..31cee41b 100644 --- a/R/class-distributionmodel.R +++ b/R/class-distributionmodel.R @@ -659,9 +659,9 @@ DistributionModel <- R6::R6Class( #' Has a offset been used? #' @return A [`logical`] flag. has_offset = function(){ - model <- self$model$offset - if(!is.Waiver(model$offset)) return( TRUE ) + if(!is.Waiver(self$model$offset)) return( TRUE ) # Also check whether offset is somehow in the equation + fit <- self$model$biodiversity[[1]]$equation ind <- attr(stats::terms.formula(fit$get_equation()), "offset") if(!is.null(ind)) return( TRUE ) }, diff --git a/R/engine_inla.R b/R/engine_inla.R index e4b866f7..11d0a46e 100644 --- a/R/engine_inla.R +++ b/R/engine_inla.R @@ -152,6 +152,7 @@ engine_inla <- function(x, mesh <- optional_mesh # Convert the study region region.poly <- methods::as(sf::st_geometry(x$background), "Spatial") + region.poly$weight <-1 # Security check for projection and if not set, use the one from background if(is.null(mesh$crs)) mesh$crs <- sp::CRS( sp::proj4string(region.poly) ) @@ -207,6 +208,7 @@ engine_inla <- function(x, # Convert the study region region.poly <- methods::as(sf::st_geometry(model$background), "Spatial") + region.poly$weight <-1 # Convert to boundary object for later suppressWarnings( @@ -219,16 +221,6 @@ engine_inla <- function(x, bdry$loc <- INLA::inla.mesh.map(bdry$loc) # Try and infer mesh parameters if not set - - # Get all coordinates of observations - locs <- collect_occurrencepoints(model, include_absences = FALSE, - tosf = FALSE) - locs <- locs[,c("x","y")]# Get only the coordinates - - assertthat::assert_that( - nrow(locs)>0, ncol(locs)==2 - ) - if(is.null(params$max.edge)){ # A good guess here is usally a max.edge of between 1/3 to 1/5 of the spatial range. max.edge <- c(diff(range(locs[,1]))/(3*5) , diff(range(locs[,1]))/(3*5) * 2) @@ -255,30 +247,42 @@ engine_inla <- function(x, params$cutoff <- cutoff } - suppressWarnings( - mesh <- INLA::inla.mesh.2d( - # Point localities - loc = locs, - # Boundary object - boundary = bdry, - # Mesh Parameters - max.edge = params$max.edge, - offset = params$offset, - cutoff = params$cutoff, - # Define the CRS - crs = bdry$crs - ) + # Get all coordinates of observations + locs <- collect_occurrencepoints(model, include_absences = FALSE, + tosf = FALSE) + locs <- locs[,c("x","y")]# Get only the coordinates + + assertthat::assert_that( + nrow(locs)>0, ncol(locs)==2 ) - # Calculate area - # ar <- suppressMessages( - # suppressWarnings( - # mesh_area(mesh = mesh, region.poly = region.poly, variant = params$area) - # ) + + # suppressWarnings( + # mesh <- INLA::inla.mesh.2d( + # # Point localities + # loc = locs, + # # Boundary object + # boundary = bdry, + # # Mesh Parameters + # max.edge = params$max.edge, + # offset = params$offset, + # cutoff = params$cutoff, + # # Define the CRS + # crs = bdry$crs + # ) # ) - # 06/01/2023: This should work and is identical to inlabru::ipoints - ar <- suppressWarnings( - Matrix::diag( INLA::inla.mesh.fem(mesh = mesh)[[1]] ) + # Convert to mesh using fmesher + mesh <- fmesher::fm_mesh_2d_inla(loc = locs, + boundary = region.poly, + offset = params$offset, + cutoff = params$cutoff, + max.edge = params$max.edge, + crs = sf::st_crs(model$background) ) + if(is.na(mesh$crs) || is.null(mesh$crs)) fmesher::fm_crs(mesh) <- sf::st_crs(x$background) + + # Calculate area weight + ar <- inlabru::fm_int(mesh)$weight |> as.vector() + assertthat::assert_that(length(ar) == mesh$n) # Now set the output diff --git a/R/engine_inlabru.R b/R/engine_inlabru.R index a482cf29..02087d01 100644 --- a/R/engine_inlabru.R +++ b/R/engine_inlabru.R @@ -137,15 +137,13 @@ engine_inlabru <- function(x, # Load a provided on mesh <- optional_mesh # Security check for projection and if not set, use the one from background - if(is.null(mesh$crs)) mesh$crs <- sp::CRS( sp::proj4string(region.poly) ) + if(is.na(mesh$crs) || is.null(mesh$crs)) fmesher::fm_crs(mesh) <- sf::st_crs(x$background) # Convert the study region region.poly <- methods::as(sf::st_geometry(x$background), "Spatial") # Calculate area - ar <- suppressWarnings( - mesh_area(mesh = mesh, region.poly = region.poly, variant = area) - ) + ar <- inlabru::fm_int(mesh)$weight |> as.vector() } else { mesh <- new_waiver() ar <- new_waiver() @@ -186,37 +184,31 @@ engine_inlabru <- function(x, # Convert the study region region.poly <- methods::as(sf::st_geometry(model$background), "Spatial") - - # Convert to boundary object for later - suppressWarnings( - bdry <- INLA::inla.sp2segment( - sp = region.poly, - join = TRUE, - crs = INLA::inla.CRS(projargs = sp::proj4string(region.poly)) - ) - ) - bdry$loc <- INLA::inla.mesh.map(bdry$loc) + region.poly$weight <-1 # Get all coordinates of observations locs <- collect_occurrencepoints(model, include_absences = FALSE, - tosf = FALSE) + tosf = TRUE) locs <- locs[,c("x","y")] # Take only the coordinates # Try and infer mesh parameters if not set if(is.null(params$max.edge)){ # A good guess here is usally a max.edge of between 1/3 to 1/5 of the spatial range. - max.edge <- c(diff(range(locs[,1]))/(3*5) , diff(range(locs[,1]))/(3*5) * 2) + max.edge <- c(diff(range(locs[["x"]]))/(3*5) , diff(range(locs[["y"]]))/(3*5) * 2) params$max.edge <- max.edge } + + # Try and infer parameter offset if(is.null(params$offset)){ # Check whether the coordinate system is longlat - if( sf::st_is_longlat(bdry$crs) ){ + if( sf::st_is_longlat(region.poly) ){ # Specify offset as 1/100 of the boundary distance - offset <- c( diff(range(bdry$loc[,1]))*0.01, - diff(range(bdry$loc[,1]))*0.01) + offset <- c( diff(range(sp::coordinates(region.poly)))*0.01, + diff(range(sp::coordinates(region.poly)))*0.01) } else { - offset <- c( diff(range(bdry$loc[,1]))*0.01, - diff(range(bdry$loc[,1]))*0.01) + warning("Default offset parameter for mesh likely won't work. Specify!") + offset <- c( diff(range(sp::coordinates(region.poly)))*0.01, + diff(range(sp::coordinates(region.poly)))*0.01) } params$offset <- offset } @@ -225,35 +217,25 @@ engine_inlabru <- function(x, # Specify as minimum distance between y coordinates # Thus capturing most points on this level # otherwise set to default - val <- min(abs(diff(locs[,2]))) + val <- min(abs(diff(locs[["y"]]))) cutoff <- ifelse(val == 0, 1e-12, val) params$cutoff <- cutoff } - suppressWarnings( - mesh <- INLA::inla.mesh.2d( - # Point localities - loc = locs, - # Boundary object - boundary = bdry, - # Mesh Parameters - max.edge = params$max.edge, - offset = params$offset, - cutoff = params$cutoff, - # Define the CRS - crs = bdry$crs - ) - ) + # Convert to mesh using fmesher + mesh <- fmesher::fm_mesh_2d_inla(loc = locs, + boundary = region.poly, + offset = params$offset, + cutoff = params$cutoff, + max.edge = params$max.edge, + crs = sf::st_crs(model$background) + ) + + # If crs is messing, use the one from background + if(is.na(mesh$crs) || is.null(mesh$crs)) fmesher::fm_crs(mesh) <- sf::st_crs(x$background) + # Calculate area - # ar <- suppressMessages( - # suppressWarnings( - # mesh_area(mesh = mesh, region.poly = region.poly, variant = params$area) - # ) - # ) - # 06/01/2023: This should work and is identical to inlabru::ipoints - ar <- suppressWarnings( - inlabru::ipoints(samplers = mesh)$weight |> as.vector() - ) + ar <- inlabru::fm_int(mesh)$weight |> as.vector() assertthat::assert_that(length(ar) == mesh$n) # Now set the output @@ -274,17 +256,16 @@ engine_inlabru <- function(x, alpha = 2, dims = c(300, 300) ) - # Convert to raster stack - out <- terra::rast( - sp::SpatialPixelsDataFrame( sp::coordinates(out), data = as.data.frame(out), - proj4string = terra::crs(self$get_data('mesh')$crs) ) - ) - terra::plot(out[[c('sd','sd.dev','edge.len')]], - col = c("#00204D","#00336F","#39486B","#575C6D","#707173","#8A8779","#A69D75","#C4B56C","#E4CF5B","#FFEA46") - ) + # Convert to raster stack + out <- stars::st_rasterize(out |> sf::st_as_sf()) |> terra::rast() + terra::plot(out,col = ibis_colours$ohsu_palette) } else { - INLA:::plot.inla.mesh( self$get_data('mesh') ) + if(inherits(self$get_data('mesh'), "fm_mesh_2d")){ + self$get_data('mesh') + } else { + INLA:::plot.inla.mesh( self$get_data('mesh') ) + } } },overwrite = TRUE) @@ -391,12 +372,10 @@ engine_inlabru <- function(x, if(mode == 'cp'){ # Create integration points by using the mesh as sampler suppressWarnings( - ips <- inlabru::ipoints( - samplers = self$get_data('mesh') - ) + ips <- inlabru::fm_int(self$get_data('mesh')) ) # Extract predictors add to integration point data - d <- get_rastervalue(coords = ips@coords, + d <- get_rastervalue(coords = sf::st_coordinates(ips), env = model$predictors_object$get_data(df = FALSE), rm.na = FALSE) for (cov in model$predictors_names) ips@data[,cov] <- d[,cov] @@ -520,6 +499,8 @@ engine_inlabru <- function(x, # include = include[[i]], # Don't need this as all variables in equation are included options = o ) + # MJ 14/70 addition. Manually conversion of likelihood data to sf + lh$data <- lh$data |> sf::st_as_sf() } # Add to list lhl[[j]] <- lh @@ -827,15 +808,14 @@ engine_inlabru <- function(x, # --- # # Fitting bru model - try({ - fit_bru <- inlabru::bru(components = comp, + fit_bru <- try({inlabru::bru(components = comp, likelihoods, options = options) }, silent = FALSE) # --- # # Security checks - if(!exists("fit_bru")){ + if(inherits(fit_bru, "try-error")){ stop('Model did not converge. Try to simplify structure and check priors!') } if(is.null(fit_bru$names.fixed)) stop('Model did not converge. Try to simplify structure and check priors!') @@ -844,20 +824,11 @@ engine_inlabru <- function(x, # Messenger if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Estimation]','green','Starting prediction.') - covs <- model$predictors_object$get_data(df = FALSE) - covs <- covs[[ which(names(covs) %in% fit_bru$names.fixed) ]] - - # Build coordinates - suppressWarnings( - preds <- inla_predpoints(mesh = self$get_data('mesh'), - background = model$background, - cov = covs, - proj_stepsize = self$get_data('proj_stepsize'), - spatial = TRUE - ) - ) + # Get Predictors as data.frame + preds <- model$predictors_object$get_data(df = TRUE,na.rm = FALSE) # Clamp? - if( settings$get("clamp") ) preds@data <- clamp_predictions(model, preds@data) + if( settings$get("clamp") ) preds <- clamp_predictions(model, preds) + preds <- sf::st_as_sf(preds, coords = c("x","y"),crs = sf::st_crs(model$background)) # Set target variables to bias_value for prediction if specified if(!is.Waiver(settings$get('bias_variable'))){ @@ -903,7 +874,10 @@ engine_inlabru <- function(x, ii <- "Intercept" } assertthat::assert_that(all( vn %in% names(preds) )) + # Set to NA all variable + preds[which(apply(preds[,vn], 1, function(z) any(is.na(z)))),vn] <- NA preds <- subset(preds, select = vn ) + # Add offset if set if(!is.Waiver(model$offset)){ ovn <- "spatial_offset" @@ -933,18 +907,21 @@ engine_inlabru <- function(x, ) pred_bru$cv <- pred_bru$sd / pred_bru$mean # Convert to raster if not set already - if(!is.Raster(pred_bru)) pred_bru <- terra::rast(pred_bru) - - # Get only the predicted variables of interest - if(utils::packageVersion("inlabru") <= '2.5.2'){ - # Older version where probs are ignored - prediction <- subset(pred_bru, c("mean","sd","q0.025", "median", "q0.975", "cv")) - names(prediction) <- c("mean","sd","q0.025", "median", "q0.975", "cv") - } else { - prediction <- subset(pred_bru, c("mean","sd","q0.05", "q0.5", "q0.95", "cv") ) - names(prediction) <- c("mean", "sd", "q05", "q50", "q95", "cv") + if(!is.Raster(pred_bru)){ + prediction <- c( + terra::rasterize(pred_bru, model_to_background(model), "mean"), + terra::rasterize(pred_bru, model_to_background(model), "sd"), + terra::rasterize(pred_bru, model_to_background(model), "q0.05"), + terra::rasterize(pred_bru, model_to_background(model), "q0.5"), + terra::rasterize(pred_bru, model_to_background(model), "q0.95"), + terra::rasterize(pred_bru, model_to_background(model), "cv") + ) + names(prediction) <- c("mean","sd","q0.05", "q0.5", "q0.95", "cv") } + # Clip to remove NA estimates + prediction <- terra::mask(prediction, model$background) + } else { prediction <- NULL } diff --git a/R/train.R b/R/train.R index 2b05f833..afeae414 100644 --- a/R/train.R +++ b/R/train.R @@ -1786,11 +1786,13 @@ methods::setMethod( if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Done]','green',paste0('Completed after ', round( as.numeric(out$settings$duration()), 2),' ',attr(out$settings$duration(),'units') )) # Quick check that the prediction is consistent with background extent - if(!is_comparable_raster(out$get_data(), x$background )){ + if(settings$get('inference_only')==FALSE){ + if(!is_comparable_raster(out$get_data(), x$background )){ o <- out$get_data() o <- terra::extend(o, x$background) |> terra::crop(x$background) out <- out$set_data("prediction", o) try({ rm(o) }) + } } # Clip to limits again to be sure diff --git a/R/utils-predictors.R b/R/utils-predictors.R index 4587d336..a350ef12 100644 --- a/R/utils-predictors.R +++ b/R/utils-predictors.R @@ -611,21 +611,40 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab is.vector(int_variables), msg = "Provide a vector of 2 variables for an interaction!") - # Make unique combinations - if(any(is.factor(env[[int_variables]]))){ - # Factor to variable interaction - i <- which(is.factor(env[[int_variables]])) - o <- explode_factorized_raster(env[[int_variables]][[i]]) - if(length(i)==length(int_variables)){ - ind <- utils::combn(names(o), 2) + # For raster or stars + if(is.Raster(env)){ + # Make unique combinations + if(any(is.factor(env[[int_variables]]))){ + # Factor to variable interaction + i <- which(is.factor(env[[int_variables]])) + o <- explode_factorized_raster(env[[int_variables]][[i]]) + if(length(i)==length(int_variables)){ + ind <- utils::combn(names(o), 2) + } else { + # Create two-way interactions as specified + ind <- rbind( matrix(names(o),nrow = 1), int_variables[which(!is.factor(env[[int_variables]]))]) + } + suppressWarnings( env <- c(env, o)) # Add to stack. } else { - # Create two-way interactions as specified - ind <- rbind( matrix(names(o),nrow = 1), int_variables[which(!is.factor(env[[int_variables]]))]) + # All possible ones + ind <- utils::combn(int_variables, 2) } - suppressWarnings( env <- c(env, o)) # Add to stack. } else { - # All possible ones - ind <- utils::combn(int_variables, 2) + if(any(Reduce("c", lapply(env_list, is.factor)))){ + i <- Reduce("c", lapply(env_list, function(z) which(terra::is.factor(z))) ) + if(length(i)>0) o <- explode_factorized_raster(env_list[[i]]) + if(length(i)==length(int_variables)){ + ind <- utils::combn(int_variables, 2) + } else { + ind <- expand.grid(matrix(unique(names(o)),nrow = 1),int_variables) |> t() + ind <- ind[,-which(apply(ind, 2, anyDuplicated)>0)] # Remove all duplicates + } + env_list[[i]] <- c(env_list[[i]], o) + rm(o) + } else { + # All possible ones + ind <- utils::combn(int_variables, 2) + } } if(is.Raster(env)){ diff --git a/tests/testthat/test_functions.R b/tests/testthat/test_functions.R index 6ff0520a..05ea1086 100644 --- a/tests/testthat/test_functions.R +++ b/tests/testthat/test_functions.R @@ -96,6 +96,15 @@ test_that('Custom functions - Test gridded transformations and ensembles', { suppressWarnings( t5 <- predictor_derivate(s, option = "int",int_variables = c("lyra", "lyrc")) ) expect_s4_class(t5, "SpatRaster") + # Make factor layer and check factor levels + ff <- terra::classify(r3, c(terra::global(r3,"min")[,1], + terra::global(r3,"mean")[,1], + terra::global(r3,"max")[,1]) ) + names(ff) <- "fact" + suppressWarnings( tff <- predictor_derivate(c(s,ff), option = "int", + int_variables = c("lyra", "fact")) ) + expect_s4_class(tff, "SpatRaster") + # --- # # Finally do some ensemble calculations ex <- ensemble(r1, r2, r3, layer = "lyr.1") diff --git a/tests/testthat/test_trainOtherEngines.R b/tests/testthat/test_trainOtherEngines.R index 76d74028..91a34759 100644 --- a/tests/testthat/test_trainOtherEngines.R +++ b/tests/testthat/test_trainOtherEngines.R @@ -518,7 +518,8 @@ test_that('Train a distribution model with INLABRU', { x <- distribution(background) |> add_biodiversity_poipo(virtual_points, field_occurrence = 'Observed', name = 'Virtual points') |> add_predictors(predictors, transform = 'none',derivates = 'none') |> - engine_inlabru() + engine_inlabru() |> + add_latent_spatial() # Train the model suppressWarnings( @@ -621,7 +622,7 @@ test_that('Train a distribution model with SCAMPR', { # Train the model suppressWarnings( - mod <- train(x, "test", inference_only = FALSE, only_linear = TRUE, + mod <- train(x |> add_latent_spatial(), "test", inference_only = FALSE, only_linear = TRUE, verbose = FALSE) )