Varying slopes LGCP model?

[Gabrielle Koerich]

Hi all,

I have been using a LGCP model for individual species distributions using inlabru. But I'm interested in now trying to fit a varying intercept and slopes model for all the species -- i.e. I want to get slopes for each covariate for each species.

However, I am not sure how this would work in an LGCP model (or even if it's possible). I believe I was successful in implementing a random intercept model, but the part that is confusing me is the random slopes, especially because I have more than one covariate (7 in fact). So far, I have tried:

How I coded the species:

occ$sp_id <- as.integer(as.factor(occ$scientific_name))
occ$sp_id1 <- as.integer(as.factor(occ$scientific_name))
N <- length(unique(occ$scientific))

The formulation:

fm <- geometry ~ Intercept(1) +
    maxT(rasters$max_temp, model = "linear") +
    degDayS(rasters$degday_season, model = "linear") +
    field(geometry, model = spde) +
    species(sp_id, model = "iid2d", n = 2*N) +
    maxT1(geometry, rasters$max_temp, copy = "species")

But I get an error:

Error in Ops.sfc(loc, knots[split]): operation >= not supported

I have also tried 

maxT1(sp_id1, rasters$max_temp, copy = "species")

Which runs, but I'm not sure this is doing what I want it to be doing (it makes more sense to me to have the response variable in this term rather than the species identity...)

So my questions are:

1. Is it possible to fit this type of model with a LGCP? or would it be preferable to fit a marked ppm?

2. If it is possible to use this formulation, how would I include slopes for more than one covariate? would it be by including + degDayS(...) + so forth?

And as a clarification, I'm not looking to estimate any covariance between species (not interested in "species interactions").

Many thanks!

[Helpdesk (Haavard Rue)]

will this tutorial help?

vignette("svc", package="INLA")



On Thu, 2024-07-11 at 02:51 -0700, 'Gabrielle Koerich' via R-inla discussion

> --
> 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 on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/00eed1ec-b5d1-4d57-9c72-0bfa8e220d3en%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Finn Lindgren]

Hi,

There are two different aspects to this:
1. spatially varying slopes
2. Species specific effects.

For 2, one approach to build the model on the product space geometry +
species, with something like
domain=list(geometry=mesh, sp_id=seq_len(number_of_species))
You didn't show enough of your code to see if you took that into
account, but the iid2d model specification looks odd; the point
process intensity needs to be a function of the domain coordinates.
But it depends on what model you want to implement. E.g., if there's a
common intensity, and you then model species _conditionally_ on the
observed point pattern, then you can have a component that depends on
sp_id,
without having to include it in the point process domain.

So the correct code depends on what model you're actually trying to
implement; it's not clear what "this type of model" actually is!
It's also essential to provide the like() calls you're using, as it's
the combination of model components and predictor expressions&model
families that define the model in the code.

For 1, I don't see a problem with the code as written; to help debug,
please run the model without the species component, to see if you
still get an error; if you do, then there's a problem with how you've
written the spatially varying coefficient definitions.
Please run traceback() after the error to help identify what part of
the code is throwing the error.
For inlabru, the vignette to look at for spatially varying coefficient
models is at https://inlabru-org.github.io/inlabru/articles/svc.html
The INLA package version of the vignette isn't helpful for inlabru
(and also hasn't been updated in a while).

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 on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/00eed1ec-b5d1-4d57-9c72-0bfa8e220d3en%40googlegroups.com.



--

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

[Gabrielle Koerich]

Hi both,

Thank you the explanation!

Essentially, I want to allow each species in my model to have their own response to each covariate, as I have observed in the individual species models that they can have different responses to the same environmental driver. So I believe that I'm not looking for a SVC where the relationship between the response variable and covariates is not the same across space -- it is not the same across species.

For 2, one approach to build the model on the product space geometry +
species, with something like
domain=list(geometry=mesh, sp_id=seq_len(number_of_species))

I have not included that; so far my model call is:

mod1 <- lgcp(fm, occ, samplers = iceFree1km,
              domain = list(geometry = mesh1),
              options = list(control.compute =
                               list(waic = TRUE,
                                    cpo = TRUE,
                                    return.marginals.predictor = TRUE),
                             control.inla = list(strategy = "simplified.laplace", int.strategy = 'eb'),
                             verbose = FALSE, safe = TRUE))

the iid2d model specification looks odd; the point
process intensity needs to be a function of the domain coordinates.

Yes, for the iid2d model specification, I based it off this thread: https://groups.google.com/g/r-inla-discussion-group/c/oVLaUgyCkwM/m/H5jR1jhJRfQJ

But I got stuck when trying to implement that for different covariates in a point intensity model.

to help debug, please run the model without the species component, to see if you
still get an error; if you do, then there's a problem with how you've written the spatially varying coefficient definitions.

I don't get an error when I run it without the maxT1(geometry, rasters$max_temp, copy = "species") component or without the species(sp_id, model = "iid2d", n = 2*N) component. So it is something with the formulation.

Complete code for a streamlined view:
spde <- inla.spde2.pcmatern(mesh = mesh1,
                              prior.range = c(100, 0.5),
                              prior.sigma = c(0.1, 0.01))

fm <- geometry ~ Intercept(1) +
    maxT(rasters$max_temp, model = "linear") +
    degDayS(rasters$degday_season, model = "linear") +
    field(geometry, model = spde) +
    species(sp_id, model = "iid2d", n = 2*N) +
    maxT1(geometry, rasters$max_temp, copy = "species")

mod1 <- lgcp(fm, occ, samplers = iceFree1km,
              domain = list(geometry = mesh1),
              options = list(control.compute =
                               list(waic = TRUE,
                                    cpo = TRUE,
                                    return.marginals.predictor = TRUE),
                             control.inla = list(strategy = "simplified.laplace", int.strategy = 'eb'),
                             verbose = FALSE, safe = TRUE))

Many thanks!

[Gabrielle Koerich]

Forgot to mention: running traceback() returns "No traceback available" but the problem is my formulation in this case.

[Finn Lindgren]

OK,

now I see more what you're trying to do. The issue is that the "copy"
feature only works if the "main" model takes the same inputs as the
"copy".
In your specification,

species(sp_id, model = "iid2d", n = 2*N) +
maxT1(geometry, rasters$max_temp, copy = "species")

where maxT1 is a spatially indexed model with spatial weights given by
the max_temp raster. But species is not a spatially indexed model, so
that won't work.
In the thread you linked they're doing something else, so it can't be
used as a blueprint here.

If all you want is to have different linear covariate effects for each
species, I believe you could do something like this:

geometry + sp_id ~ 0 +
Species_Intercept(sp_id, model="iid") +
maxT(rasters$max_temp, model="linear", group=sp_id,
control.group=list(model="exchangeable")) +
degDay(...) +
...

and so on for each spatial covariate. Since your model will need to
have point pattern intensity lambda(space, species) (different
intensity for each species), you'll need
domain=list(geometry=mesh, sp_id=seq_len(N))
in the lgcp() call.

The above approach can be combined with having a "common effect" of
each covariate as well, but for spatially constant effects that
probably wouldn't make much of a difference to the estimates.

Add
+ field(geometry, model=spde)
for a common random field effect, or do the same with that one;
+ field(geometry, model=spde, group=sp_id, ...)

Optionally, you could try
+ common_field(geometry, model=spde_common) +
species_field(geometry, model=spde_species, group=sp_id, ...)
where optionally you'd have used constr=TRUE when creating
spde_species, to force an integrate-to-zero constraint, which might
help with identifiability against the common_field component (but only
globally; they wouldn't be identifiable locally, so might have to be
very specific with how you set the priors for the correlation range to
nudge the components to capture what you want them to be able to
capture.

Finn

On Thu, 11 Jul 2024 at 13:39, 'Gabrielle Koerich' via R-inla

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e79da9a6-61cb-47af-96a4-e555c9780968n%40googlegroups.com.

[Gabrielle Koerich]

This is exactly what I was looking for, and it makes way more sense! Thank you so much.

I'm running some tests and will try out including a common random field.

Cheers!

[Gabrielle Koerich]

Hi again,

I have been trying out the new model specification and I have managed to fit models with only the SPDE term, after adjusting some of the priors (the iid and fixed effects) and adding constr = T to the SPDE:

spde <- inla.spde2.pcmatern(mesh = mesh,
                                prior.range = c(130, 0.01),
                                prior.sigma = c(0.01, 0.001),
                                constr = T)

fm <- geometry + sp_id ~ 0 +
  speciesIntercept(sp_id, model = "iid", hyper=list(prec=list(prior="pc.prec", param=c(0.1, 0.001)))) +
  field(geometry, model = spde)

mod <- lgcp(fm, occ, samplers = iceFree1km,
            domain = list(geometry = mesh, sp_id = seq_len(N)),

            options = list(control.compute =
                             list(waic = TRUE,
                                  cpo = TRUE,
                                  return.marginals.predictor = TRUE),
                           control.inla = list(strategy = "simplified.laplace", int.strategy = 'eb'),

                           control.fixed = list(prec.intercept = 1, prec = 1),
                           verbose = F, safe = TRUE))

summary(mod)
inlabru version: 2.11.1
INLA version: 24.06.27
Components:
speciesIntercept: main = iid(sp_id), group = exchangeable(1L), replicate = iid(1L)
field: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
Likelihoods:
  Family: 'cp'
    Data class: 'sf', 'data.frame'
    Predictor: geometry + sp_id ~ .
Time used:
    Pre = 0.766, Running = 76, Post = 5.35, Total = 82.2
Random effects:
  Name   Model
    speciesIntercept IID model
   field SPDE2 model

Model hyperparameters:
                                  mean    sd 0.025quant 0.5quant 0.975quant    mode
Precision for speciesIntercept   0.171 0.015      0.143    0.170      0.202   0.169
Range for field                140.095 6.865    127.056  139.931    154.078 139.611
Stdev for field                  0.926 0.026      0.876    0.926      0.978   0.925

Deviance Information Criterion (DIC) ...............: -51145.53
Deviance Information Criterion (DIC, saturated) ....: -51693.79
Effective number of parameters .....................: -56347.89

Watanabe-Akaike information criterion (WAIC) ...: 5367.01
Effective number of parameters .................: 86.59

Marginal log-Likelihood:  -32381.66
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

But when I try adding the covariates (with the same SPDE):

fm1 <- geometry + sp_id ~ 0 +
      speciesIntercept(sp_id, model = "iid", hyper=list(prec=list(prior="pc.prec", param=c(0.1, 0.001)))) +
      twi(rasters$twi, model = "linear", group = sp_id, control.group=list(model="exchangeable")) +
      field(geometry, model = spde)

mod1 <- lgcp(fm, occ, samplers = iceFree1km,

                  domain = list(geometry = mesh, sp_id = seq_len(N)),

                  options = list(control.compute =
                                   list(waic = TRUE,
                                        cpo = TRUE,
                                        return.marginals.predictor = TRUE),
                                 control.inla = list(strategy = "simplified.laplace", int.strategy = 'eb'),

                                 control.fixed = list(prec.intercept = 1, prec = 1),
                                 verbose = T, safe = TRUE))

I keep getting this error:

*** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1
Error: C stack usage  64282397 is too close to the limit

It shouldn't be a memory issue as it looks like the model is using about 50% of the memory in my computer (which has 32gb), it seems to be a convergence issue. I'd like to know if there are any suggestions on this or if it's just a problem with the data/complexity of the model. I have also tried adding a constr = T to the iid term (speciesIntercept) to see if it would change anything, but that didn't help.

I have also added the mod$logfile here for the model without covariates for reference. I am running it on R 4.4.1 and the INLA testing version.

Many thanks!

[Finn Lindgren]

Hi,

I don’t see anything immediately wrong. Any chance you can capture the verbose on console output for the run that fails? Hard to debug without info.

Finn

On 18 Jul 2024, at 16:44, 'Gabrielle Koerich' via R-inla discussion group <r-inla-disc...@googlegroups.com> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/936e7974-9341-4744-940a-e1f59dfc237fn%40googlegroups.com.

<mod.txt>

[Gabrielle Koerich]

Hi Finn,

Unfortunately not, I get only the error message. Traceback shows this:

> traceback()
6: .makeMessage(..., domain = domain, appendLF = FALSE)
5: data.frame(message = .makeMessage(..., domain = domain, appendLF = FALSE),
       timestamp = Sys.time(), verbosity = as.integer(verbosity))
4: bru_log_message(paste0("iinla: Problem in inla: ", result), verbose = FALSE,
       verbose_store = options$bru_verbose_store)
3: iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]])
2: bru(components, lik, options = options, .envir = .envir)
1: lgcp(fm1, occ, samplers = iceFree1km, domain = list(geometry = mesh,

       sp_id = seq_len(N)), options = list(control.compute = list(waic = TRUE,
       cpo = TRUE, return.marginals.predictor = TRUE), control.inla = list(strategy = "simplified.laplace",
       int.strategy = "eb"), control.fixed = list(prec.intercept = 1,
       prec = 1), verbose = T, safe = TRUE))

[Finn Lindgren]

That’s strange; no output at all except for the error message, event thought you set verbose=TRUE?

Finn

On 18 Jul 2024, at 19:01, 'Gabrielle Koerich' via R-inla discussion group <r-inla-disc...@googlegroups.com> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/c60ffa39-4a7f-4398-af35-58c8b70aa3ffn%40googlegroups.com.

[Finn Lindgren]

The “stack error” is usually a side-effect of the error reporting, and not the actual problem.

Finn

On 18 Jul 2024, at 19:21, Finn Lindgren <finn.l...@gmail.com> wrote:

That’s strange; no output at all except for the error message, event thought you set verbose=TRUE?

[Gabrielle Koerich]

Hi Finn,

Thank you very much for the direction! I tried to get a message error by reducing the number of species to 2 and figured out what it's likely happening.

I got a warning saying that there were NAs in the covariates and that inlabru was trying to fill those NA values. So probably the "stack error" showed up before any verbose from the model because the bru_eval_spatial would get stuck before the model initialisation. I think the sf package updated how it deals with crs reprojections recently because I had some errors with the way I was reprojecting spatial objects (and that's what gave me the clue to rerun the mesh code which led me to the warning). I'll fix that and I believe that will solve the issue. Thanks!

[Finn Lindgren]

Hi,

I wasn’t aware of any changes in reprojections, so would be keen to know more about that issue..

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/86875c31-b30a-468e-a69b-ab02d701af61n%40googlegroups.com.

[Gabrielle Koerich]

I'm back -- I made sure everything is correctly in sf format and also that the covariates are no longer being filled by inlabru, and I still get the same "stack error" message without getting an output before it crashes even though I am setting verbose = T.

Running a model with the covariate but no species-specific term in the covariate works:

fm1 <- geometry + sp_id ~ 0 +

      speciesIntercept(sp_id, model = "iid") +
      twi(rasters$twi, model = "linear") +
      field(geometry, model = spde)

mod1 <- lgcp(fm1, occ, samplers = iceFree1km,

                  domain = list(geometry = mesh, sp_id = seq_len(N)),
                  options = list(control.compute =
                                   list(waic = TRUE,
                                        cpo = TRUE,
                                        return.marginals.predictor = TRUE),
                                 control.inla = list(strategy = "simplified.laplace", int.strategy = 'eb'),

                                 verbose = T, safe = T))

And I get a model that looks ok (verbose/summary output attached), and confirms that my covariates are complete now.

So it seems like the problem is when I add the species-specific term, unfortunately:

twi(rasters$twi, model = "linear", group = sp_id, control.group=list(model="exchangeable"))

I'm sorry for not being able to get a more useful error output (even reducing it to two species doesn't work).

[Finn Lindgren]

Hi,
problem identified; I wasn't aware that model="linear" isn't allowed
with 'group'.
I get the following warning when running this simple construction:

> fit<-bru(y~0+eff(covar,model="linear",group=id),data=data.frame(y=rnorm(4),covar=1:4,id=c(1,2,1,2)))

*** inla.core.safe: The inla program failed, but will rerun in case
better initial values may help. try=1/1

Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]],
options = info[["options"]]) :
iinla: Problem in inla: Error in inla.core.safe(formula = formula,
family = family, contrasts = contrasts, :
Model = 'linear' do not accept argument 'group'. Please use another
model component.


So a solution is to convert model="linear" into a different model, but
set it up to be effectively the same thing.
Here's one way to do it (not sure of the default precision for "fixed"
effects). It's probably better to turn the
"group" effect into the "main" effect instead, but this was the most
similar alternative I could come up with on short notice...

> fit<-bru(y~0+eff(1,weights=covar,model="iid",group=id,
hyper=list(prec=list(initial=log(1/30^2), fixed=TRUE))),
data=data.frame(y=rnorm(4),covar=1:4,id=c(1,2,1,2)))

Finn

On Fri, 19 Jul 2024 at 13:29, 'Gabrielle Koerich' via R-inla

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/6d6684d7-8344-4071-a0eb-2b7872d3a3f1n%40googlegroups.com.

[Gabrielle Koerich]

Oh I see! This works, thank you very much. I'll try to understand better what this is doing but it seems to work fine for now.

Cheers!