Help with Data Integration and an issue related to the background layer

[Daviser95]

Hi!
I am recently using some R packages to perform SDMs, and in the last days, I have been trying the ibis.iSDM.
I think this package would be the best option for what I want to do, but I am having some trouble for integrating data.
Specifically, I have a large GPS-based with ~90.000 occurrences for a large study area. I decided to split the dataset into 10 smaller datasets, and I want to use nine to build models and one to validate them.

The first option I considered was to build different models and then combine them all to create an ensemble.
But after reading the paper and the package functionalities, I would like to create an iterative process in which the output from the first modeling process with the 'dataset1', is used to improve the prediction from the modeling procedure with the 'dataset2'. Then, this new output is used for the modeling procedure with 'dataset3' and so on.

However, I am not sure how to integrate this data. Currently, I am using the option 'add_predictors_model', but this uses the previously obtained model as a variable during the procedure, and I am not sure if this is the best option to integrate this data. I wonder if integrating this data as a prior is better, but I am not able to do this and I found only a few examples in the reference site.

Considering my objective, I would like to know if this solution could work or if there are other options that I should consider in order to iteratively integrate the models, dataset by dataset. Your opinion would be very important here.

This is my code:

library(ibis.iSDM)
library(inlabru)
library(INLA)
library(uuid)
library(assertthat)
library(Ckmeans.1d.dp)
library(pdp)
library(terra)
library(sf)
library(spThin)

background <- rast("G:/davide/variables_100m/all_final_variables/dem_100m_wgs.tif")
background_v <- st_read("G:/davide/area_study.shp")

path_vars<-"G:/davide/variables_100m/selected_current" ### qui metti il percorso alla cartella dove stanno le variabili
env <- list.files(path_vars,pattern=".tif$")
#setwd()
EnvVars<-rast(lapply(env,function(x){
  var<-rast(paste0(path_vars,"/",x))
  return(var)}))
EnvVars
names(EnvVars) <- c("bio15", "bio17", "bio18", "bio3", "bio9", 
                    "elevation", "EucDist_rivers_bigger", "dist_roads", "dist_urbanareas", "gfcc_treecov",
                    "hfp2019", "Percent_NonVegetated", "NDVI", "Npp", "slope")
#predictors <- subset(EnvVars, c("elevation", "Percent_NonVegetated", "NDVI", "slope"))

data_df <- read.csv("G:/davide/gps_data_reduced_100m.csv")

data_full <- read.csv("G:/davide/gps_data_reduced_100m.csv") %>% 
  st_as_sf(coords=c("Longitude","Latitude"), crs=4326)

shuffled_data_df= data_full[sample(1:nrow(data_df)), ]
dat1 <- shuffled_data_df[1:8750,]
dat2 <- shuffled_data_df[8751:17500,]
dat3 <- shuffled_data_df[17501:26250,]
dat4 <- shuffled_data_df[26251:35000,]
dat5 <- shuffled_data_df[35001:43750,]
dat6 <- shuffled_data_df[43751:52500,]
dat7 <- shuffled_data_df[52501:61250,]
dat8 <- shuffled_data_df[61251:70000,]
dat9 <- shuffled_data_df[70001:78750,]
  
dat_validation <- shuffled_data_df[78751:87574,]  

dat1$pres_abs <- 1
# pseudoabsences

densities_ranges <- sf::st_read("G:/davide/highdens_sites.shp")

abs <- pseudoabs_settings(background = background,
                          nrpoints = 8750, # Sample 1000 points
                          method = "range", # Option for minimum convex polygon
                          inside =  FALSE, # Sample exclusively outside
                          layer = densities_ranges,
)
print( abs )

fulldat1 <- add_pseudoabsence(dat1,
                            # Point to the column with the presence information
                            field_occurrence = 'pres_abs',
                            settings = abs)

plot(fulldat1['pres_abs'])

# predictor transformations
envar_trans <- predictor_transform(EnvVars, option = "norm") 

# Harmonize the predictors
envar_ready <- predictor_homogenize_na(env = envar_trans)

# set first model

mod_data1 <- distribution(background_v) |> 
  add_biodiversity_poipa(fulldat1,
                         name = "dataset1_2",
                         field_occurrence = "pres_abs") |>  
  add_predictors(env = envar_ready) |> 
  engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)


# Finally train
fit1 <- train(mod_data1,
             runname =  "Test xgboost run1",
             verbose = FALSE # Don't be chatty
)

# Visualize Results

# Plot the mean of the posterior predictions
plot(fit1, "mean")
writeRaster(fit1$fits$prediction, "G:/davide/Namibia/result/ibis/curr_pred_xgboost.tif")

effects(fit1)

# To calculate a partial effect for a given variable
parteff_elev <- ibis.iSDM::partial(fit1, x.var = "elevation", plot = TRUE)
# Similarly the partial effect can be visualized spatially as 'spartial'
parteff_elev_plot<- spartial(fit1, x.var = "elevation")
plot(parteff_elev_plot[[1]], col = rainbow(10), main = "Marginal effect of NDVI on the relative reporting rate")

# Calculate a threshold based on a 50% percentile criterion
fit1_threshold <- threshold(fit1, method = "percentile", value = 0.5)

# Notice that this is now indicated in the fit object
print(fit1)
# There is also a convenient plotting function
fit1_threshold$plot_threshold()

# It is also possible to use truncated thresholds, which removes non-suitable areas
# while retaining those that are suitable. These are then normalized to a range of [0-1]
fit1_threshold_trunc <- threshold(fit1, method = "percentile", value = 0.5, format = "normalize")
fit1_threshold_trunc$plot_threshold()

# Validation of model predictions

# By Default validation statistics are continuous and evaluate the predicted estimates against the number of records per grid cell.
fit1$rm_threshold()
validate(fit1, method = "cont")

# If the prediction is first thresholded, we can calculate discrete validation estimates (binary being default)
fit1_threshold_bin <- threshold(fit1, method = "percentile", value = 0.5, format = "binary")
validate(fit1_threshold_bin, method = "disc")

# repeat for Dataset2

dat2$pres_abs <- 1
sdm1 <- rast("G:/davide/result/ibis/curr_pred_xgboost.tif")

fulldat2 <- add_pseudoabsence(dat2,
                              # Point to the column with the presence information
                              field_occurrence = 'pres_abs',
                              settings = abs)
plot(fulldat2['pres_abs'])

# set model and add previous model for dat1

mod_data2 <- distribution(background_v) |> 
  add_biodiversity_poipa(fulldat2,
                         name = "dataset2",
                         field_occurrence = "pres_abs") |>  
  add_predictors(env = envar_ready) |> 
  engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)|>
  add_predictors_model(model = fit1)

#mod_2_add <- mod_data2 |>
#  add_predictors_model(model = fit1)
#p1 <- XGBPrior(variable = , hyper = "increasing", ...)
# Single/Multiple priors need to be passed to `priors` and then added to the model object.
#pp <- priors(p1)
# The variables and values in this object can be queried as well
#pp$varnames()
# Priors can then be added via 


# Finally train
fit2 <- train(mod_data2, runname =  "Test xgboost run2 plus 1",
              verbose = FALSE # Don't be chatty
)

# Visualize Results

# Plot the mean of the posterior predictions
plot(fit2, "mean")

Finally, I am getting an error when I use as a background 'background', that is a raster layer (specifically, is the elevation layer listed in the predictors). When using this layer, I get an error I can't understand:

# will give error
mod_data1_err <- distribution(background) |> 
 add_biodiversity_poipa(fulldat1,
                        name = "dataset1_2",
                        field_occurrence = "pres_abs") |>  
 add_predictors(env = envar_ready) |> 
 engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)

Errore in h(simpleError(msg, call)) : 
 errore durante la valutazione dell'argomento 'x' nella selezione di un metodo per la funzione 'add_predictors': errore durante la valutazione dell'argomento 'x' nella selezione di un metodo per la funzione 'add_biodiversity_poipa': Values in background layer can only be unique!

Thank you in advance for your help!
Best regards
Davide

[Martin-Jung]

Hi Davide,

Nice to hear that you plan to use the package!
Some quick responses:

• Generally it is hard for me to understand your issues. If you find a bug, then please create one to reproduce the issue for people without access to your input data (see http://adv-r.had.co.nz/Reproducibility.html ).
• You seem to use the engine_xgboost which only supports offsets, predictors and priors as types of integration. So you could think of either first fitting a model as you do and then set it as predictor, do a sequential prediction (adding two biodiversity datasets) or an ensemble, etc. Many different ways that can or cannot work well.
• If you want to explore some more examples using the package, there is - besides the online vignette articles - also the supplementary materials from the published paper (can be found here).

With regards to your error at the end:

• The error message at the end gives you a hint: Values in background layer can only be unique. So the issue is that your background layer (background) has more than 1 unique values, which is not allowed. Simply recategorize all values to 1 for example (background[background>0] <- 1).

As this is not really an issue with regards to the package but data preparation, I have converted this into a discussion entry (see github tab).

Cheers,
Martin

[Martin-Jung]

Hi Davide,

1. The suitability prediction output is indeed weird and it might that be that there is some bug in the projection code for engine_xgboost pending that everything else works. Could you try out
   fit8a$project(envar_mcp)
   if that works?
   Also consider installing the current dev built when in doubt (remotes::install_github("iiasa/ibis.iSDM", "dev")).
   That being said the current dev built works for me, so not quite sure what happened (without test data to reproduce the issue).
   This works:

library(ibis.iSDM)
background <- terra::rast(system.file('extdata/europegrid_50km.tif', package='ibis.iSDM'))
virtual_range <- sf::st_read(system.file('extdata/input_data.gpkg', package='ibis.iSDM'), 'range',quiet = TRUE)
virtual_points <- sf::st_read(system.file('extdata/input_data.gpkg', package='ibis.iSDM'), 'points',quiet = TRUE)
ll <- list.files(system.file('extdata/predictors_presfuture/',package = "ibis.iSDM",mustWork = TRUE),
                 full.names = T)
# Load the same files future ones
suppressWarnings(
  pred_future <- stars::read_stars(ll) |> stars:::slice.stars('Time', seq(1, 86, by = 10))
)
sf::st_crs(pred_future) <- sf::st_crs(4326)

pred_current <- ibis.iSDM:::stars_to_raster(pred_future, 1)[[1]]
names(pred_current) <- names(pred_future)

x <- distribution(background) %>%
  add_biodiversity_poipa(virtual_points |> add_pseudoabsence(field_occurrence = 'Observed',template = background),
                         field_occurrence = 'Observed', name = 'Virtual points') %>%
'binomial') %>%
  add_predictors(pred_current, transform = 'scale',derivates = "none") %>%
  engine_xgboost(iter = 5000)

modf <- train(x, runname = 'Null',only_linear = F,
              inference_only = F, verbose = T)
              
mod <- scenario(modf, limits = subs) |> 
  add_predictors(env = pred_future, transform = 'scale', derivates = "none") |>
  project()
plot(mod)

2. The scenario projections require a terra SpatRaster or a stars object as predictor. Regardless of type, the predictor (from for example WorldClim) needs to have a time dimension in Date or PosixCT format. If you don't need the date information (for example because you are predicting on a single time slot), just set it to some dummy value, e.g. 2020-01-01.

[Daviser95]

Hi Martin,

Thank you for your quick answers. I tried with fit8a$project(envar_mcp) and fit8a$project(envars_ready) but it didn't work, returning the error Error: newdata is not a data frame or newdata does not inherit from class xgb.DMatrix. I also tried installing the current dev, but nothing changed, and unfortunately now when I run:

sc1 <- scenario(fit8a) |> 
 # Apply the same variable transformations as above. 
 add_predictors(envar_ready, transform = 'scale')

It returns me the following error: Error in model$predictors object$get transformed params() : attempt to apply non-function

I don't know how could I solve these issues. I could also share the data used if you want to reproduce this issue.

Thank you!
Davide

[Martin-Jung]

Yes, the current dev built contains a fix to the predictor scaling (see issue #113 ). You would need to refit your model in this case (likely the reason why the get_transformed_params() value does not exist).

fit8a$project(envars_ready). Ahh srry, this might not accept a SpatRaster object as input but expects a data.frame. Try convertig first with terra::as.data.frame(envars_ready, xy =TRUE)

[Daviser95]

Thank you for these suggestions!
So, I made envar_ready_df <-terra::as.data.frame(envars_ready, xy =TRUE) and then try_from_Martin <- fit8a$project(envars_ready) and works, but even in this case when I plot it, the projection is biased as the image attached in the previous days.
Moreover, I transformed the predictors as follows:

envar_mcp <- predictor_transform(envar_mcp, "scale") #predictors for the calibration area (buffered mcp)
envar_ready <- predictor_transform(envar_ready, "scale") #predictors for the entire study area

then:

mod_8a <- distribution(background_mcp) |>
  add_biodiversity_poipa(fulldat_gbif, field_occurrence = "pres_abs", docheck = F) |>
  add_predictor_range(layer = densities_ranges, method = "distance", distance_max = 5e5) |>
  add_biodiversity_poipa(fulldat1_range, field_occurrence = "pres_abs",docheck = F) |>
  add_predictors(envar_mcp) |>
  engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)
fit8a2 <- train(mod_8a, method_integration = "prior", only_linear = FALSE, verbose = TRUE)

But when I try to refit the model, I get the error Error in method[i] : object of type 'environment' is not subsettable.
I don't know what is not working, in another R Project I am able to run models, for example

x1 <- distribution(bg) |>
  add_predictor_range(layer = iucn_extant_range, distance_max = 5e5) |>
  add_biodiversity_poipa(fulldata, field_occurrence = "pres_abs", docheck = F) |>
  add_predictors(predictors, transform = "scale") |>
  engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)

fit_x1 <- train(x1, only_linear = FALSE, verbose = TRUE

with no issues.
Again, thank you for your support!

[Martin-Jung]

Hi,
you are using xgboost and the method_integration = "prior" with 2 datasets. What this means is that the coefficients from the fitted first dataset are used as monotonic constraints (e.g. fix the coefficients) for the second model and dataset.
I have not encountered the error that you report before and it might be a bug. Could you try and recreate it with testing data in the package or upload me a small subset of the data you use?
(my mail is in the package DESCRIPTION file).

[Martin-Jung]

Hi,
thanks for sending me a script and testing data. The script you send me mostly works (with exception of the corruped densities shapefile and file directories) and I could reproduce the stripped prediction issue with xgboost. This has now been fixed in commit f79b88d on branch dev, so update to the dev branch to test...

Generally:

• You cannot project with a predictor set that does not include all predictor names used in the model. In your case you are using the range as well as the prediction from one dataset as input into another. Check with summary(fit) that your covariates match.
• To add (previous) prediction from another model there are some hacky things you can do for example via
  my_prediction <- fit8a$model$predictors_object$get_data()[["poipa_Species.._mean"]] but also have a look at add_predictors_model().
• Other than memory efficiency still not sure why you using one predictor set to then predict on a larger predictorset. You could avoid yourself a lot of pain by just using the large one throughout. For projection the whole raster is read into memory anyway, so you not really saving anything.

Hope that helps and good luck with your modelling exercise. Closing this now.

Appendix:

#modelling
mod_8a <- distribution(background_mcp) |>
  add_biodiversity_poipa(fullgbif, field_occurrence = "pres_abs", docheck = F) |>
  add_predictor_range(layer = densities_ranges, method = "distance", distance_max = 5e5) |>
  add_biodiversity_poipa(fulldat1, field_occurrence = "pres_abs",docheck = F) |>
  add_predictors(envar_mcp) |>
  engine_xgboost(booster = "gbtree", iter = 8000L, learning_rate = 0.001, gamma = 6)
fit8a <- train(mod_8a, method_integration = "predictor", only_linear = FALSE, verbose = TRUE)

plot(fit8a, "mean")
#aggiungo anno a variabili di scenario

# Check coefficients that are needed for any projection
summary(fit8a) # poipa_Species.._mean missing!
# This is the prediction from your previous model
my_prediction <- fit8a$model$predictors_object$get_data()[["poipa_Species.._mean"]]
# It however differs from the scale you want to project on
# Thus we need to extend it.
my_prediction <- terra::extend(my_prediction, envar_ready)

# Get "pres_abs" prediction and add to model object as otherwise missing
envar_ready <- c(envar_ready, my_prediction)
time_vect <- rep("2020-01-01", terra::nlyr(envar_ready)) |> as.Date()
terra::time(envar_ready) <- time_vect
has.time(envar_ready)

sc1 <- scenario(fit8a) |>
  # Apply the same variable transformations as above.
  add_predictors(env = envar_ready, transform = 'none')

sc.fit1 <- sc1 |> project()
sc.fit1