biodiverse / spOccupancy

Single-species, Multi-species, and Integrated Spatial Occupancy Models
https://www.jeffdoser.com/files/spoccupancy-web/
GNU General Public License v3.0
52 stars 8 forks source link

Issue with specifying detection and occupancy formulas in function iterating over a list of covariate names #27

Closed binturong123 closed 1 year ago

binturong123 commented 1 year ago

Hello, thanks for this great package. I am trying to run a list of univariate multispecies models as an array job but am having trouble given that the formula has to be specified with the exact name of the variable. Thus I am not able to iterate over my list and receive an error reporting that Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric. I've tried the "cat()" and "str2lang()" but both did not work. I would appreciate your assistance if there is any work-around with this. Cheers.

my_functiona <- function(x) {
occ.ms.sp.formula <- ~ 1 det.ms.sp.formula <- ~ scale(x)

out.sp.ms <- msPGOcc(occ.formula = occ.ms.sp.formula, det.formula = det.ms.sp.formula, data = mainframe, inits = ms.inits, priors = ms.priors, verbose =TRUE, n.samples = 50000, k.fold = 10, n.thin = 50, n.burn = 25000, n.chains = 3, n.report = 1000)

return(out.sp.ms)
}

I

doserjef commented 1 year ago

Hello!

Great question. Take a look at the example code below with the built-in hbef2015 data set in spOccupancy, which you should be able to adapt to your data set. You can use the paste() and formula() functions to accomplish what you want.

library(spOccupancy)
# Example data set
data(hbef2015)

# Specify each element of list as a set of covariates to include in a given model run
det.vars <- list(a = c('scale(day)', 'scale(tod)'), 
                 b = c('scale(tod)')) 
# List to hold model results
out.models <- vector(mode = 'list', length = length(det.vars))
for (i in 1:length(det.vars)) {
  # Paste variables together if multiple, paste with ~, then convert to formula
  curr.det.formula <- formula(paste('~', paste(det.vars[[i]], sep = '', collapse = '+')))
  # Run the model
  out.models[[i]] <- msPGOcc(occ.formula = ~ 1, 
                             det.formula = curr.det.formula, 
                             data = hbef2015, 
                             n.samples = 100, 
                             n.report = 25,
                             n.burn = 0, 
                             n.thin = 1, 
                             n.chains = 1)
}

Let me know if that doesn't solve your problem.

Jeff

binturong123 commented 1 year ago

Hi Jeff,

Amazing, thank you so much! This worked perfectly. The only thing with array jobs that may be an issue is the k.fold specs that automatically implements foreach() or mclapply(). Other than that, its working great. Thank you so much again!

Best, Arata