Problem with Spatio-Temporal Log Gaussian Cox Process Model

[Greta Panunzi]

Hello everyone, I am trying to estimate a Log Gaussian Cox Process Model in which the effort function changes every year. 

The information on effort is contained in a list of dataframes defined as follows:

List of 24
 $ eff_2000:'data.frame': 11589919 obs. of  4 variables:
  ..$ x      : num [1:11589919] -1760 -1759 -1758 -1757 -1756 ...
  ..$ y      : num [1:11589919] 5473 5473 5473 5473 5473 ...
  ..$ eff2000: num [1:11589919] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ Year   : int [1:11589919] 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
 $ eff_2001:'data.frame': 11589919 obs. of  4 variables:
  ..$ x      : num [1:11589919] -1760 -1759 -1758 -1757 -1756 ...
  ..$ y      : num [1:11589919] 5473 5473 5473 5473 5473 ...
  ..$ eff2001: num [1:11589919] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ Year   : int [1:11589919] 2001 2001 2001 2001 2001 2001 2001 2001 2001 2001 etc...

Meanwhile, the effort function is defined like this:

log_detect_year = function(xy, t1, t2, Year1) {

indx <- Year1 

 print(paste("Detection:", indx), sep = " ") 

 print(xy) 

 df <- eff.r.ls[[indx]] 

 eff <- rast(df[,1:3], type = "xyz") 

 dens = eval_spatial(eff, xy) 

 pnorm(dens / qexppnorm(t1, rate = 0.5) + t2, log.p = TRUE) 

} 

My INLA model is defined as follows:

cmp = ~

  intercept1(1, mean.linear = 0, prec.linear = 2) +
  SPDE(geometry, model = barrier.model,
       mapper = bru_mapper(mesh_medit),
       replicate = Year1, nrep = 4) +
  depth_s(depth.r$depth, mean.linear = 0.5, prec.linear = 2) +
  temp(temp.mean.r$temp_fa, mean.linear = 0, prec.linear = 1.5) +
  bpi(bpi.r$bpi, mean.linear = 0, prec.linear = 1.5) +
  Year_cov(Year1, model = "linear", prec.linear = 0.01) +
  effort_scale(1, prec.linear = 1) +
  effort_shift(1, prec.linear = 10)

formula_t = geometry + Year1 ~
  intercept1 +
  SPDE +
  depth_s +
  temp +
  bpi +
  Year_cov +
  log_detect_year(geometry, effort_scale, effort_shift, Year1)

lik_t <- like("cp",
              formula = formula_t,
              samplers = boundary_domain,
              domain = list(geometry = mesh_medit,
                            Year1 = seq(21, 24,1)),
              data = mk.sf)

fit <- bru(components = cmp,
           lik_t, 
           options = list(
             verbose = FALSE,
             inla.mode = "experimental",
             bru_max_iter = 35,
             control.inla = list(int.strategy = "eb")
           ))

The model doesn’t run and returns this error:

Error in eff.r.ls[[indx]] : no such index at level 2

I’m definitely doing something wrong in the implementation. 

Can you help me figure it out?

[Finn Lindgren]

Hi Greta,

I think there may be several issues with the log_detect_year() function, but can't easily say what the main issue is.
Also, at least the SPDE component definition doesn't match the domain definition for Year1; Either Year1 needs to be in {1,2,3,4} everywhere, or, if you wish to keep it as {21,22,23,24}, you need to specify mappers to map those values down to 1:4-indices; the replicate information needs to have indices in {1, ..., nrep} unless you supply a mapper to care of that mapping.

As Year1 takes multiple different values, you need to write log_detect_year() year in a vectorised way,

e.g., by finding out which elements should extract values for a given year;

values <- numeric(NROW(xy))

for (y in unique(Year1)) {

  in_y <- Year1 == y

  df <- eff.r.ls[[y]]

  eff <- ...

  dens <- eval_spatial(eff, xy[in_y])

  values[in_y] <- ...

}

What kind of effort is it? Could you perhaps pre-convert it to a raster with one layer per year?
Then eval_spatial() would be able to handle the vectorisation over the different years automatically (if terra has been updated with the bugfix for numerical indexing I submitted earlier this year; I haven't checked if that has made it to CRAN yet).
But if the effort is line transect effort, you could pre-convert it to an sf object with lines and linestrings perhaps.

Finn

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/529986e3-bfca-4bc6-8035-d519f97c14b4n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com