Spatio-temporal separable and non-separable LGCP with temporal variying spatial covariates

[H. Andrés Ramírez Maguiña]

Hello,

I know that similar topics have already been discussed on this forum, but I would be extremely grateful if you could guide me regarding my particular case.

I am working with epidemiological data from agricultural surveillance, where I am modeling the spread of a crop disease across a regional agricultural area. Specifically, I am trying to fit a non-separable spatio-temporal log-Gaussian Cox Process model using the diffusion extension of Gaussian Matérn Fields proposed by Lindgren et al. (2024) (https://hdl.handle.net/2117/421369), implemented through the INLAspacetime package

My response variable consists of point occurrences of the disease, including spatial coordinates and collection periods. The temporal unit corresponds to four-month periods over a 5-year time span.

I also have continuous spatial covariates that are assumed to be time-invariant, such as terrain topography, slope, among others. I incorporated them into the model following the example by Borchers and Lindgren (https://inlabru-org.github.io/inlabru/articles/2d_lgcp_covars.html), using custom functions:

formula <- geometry + t ~ Intercept(1) +
Aspect(f.Aspect(.data.), model = "linear") +
Curvate(`f.Cross-Sectional Curvature`(.data.), model = "linear") +
DEM(f.DEM(.data.), model = "linear") +
Eatness(f.Eatness(.data.), model = "linear") +
Slope(f.Slope(.data.), model = "linear") +
time = t),
model = stmodel)

For example, my custom function for curvature is defined as:

f.Curvature <- function(where) {
# Extract the values
v <- eval_spatial(Curvature, where, layer = "elevation")
# Fill in missing values; this example would work for
# SpatialPixelsDataFrame data
# if (anyNA(v)) {
# v <- bru_fill_missing(elev, where, v)
#
v
}

My covariates are vectorized rasters clipped to match the mesh area and stored in a list, for example raster_mask$Curvature.

My spatio-temporal model is defined as follows, where smesh is the spatial mesh, tmesh is the temporal mesh, and model <- "202" (non-separable) or model <- "220" (separable) . stModel.define is a function from the INLAspacetime package:

stmodel2 <- stModel.define(
smesh, tmesh, model,
control.priors = list(
prs = c(5, 0.05),
prt = c(3, 0.5),
psigma = c(2, 0.05)),
constr = TRUE

Up to this point, everything works well. However, my main question is the following: how should I incorporate spatial covariates that vary over time?  Another question would be whether the implementation would be the same for both the non-separable and the separable model.

For example, I have NDVI, temperature, and precipitation rasters corresponding to the disease sampling periods (every four months), using the same discrete temporal units as my response variable.  My dataset includes 10 sampling periods (10 four-month intervals), so for each covariate I have 10 different raster layers, one for each time period.

I would greatly appreciate any guidance, suggestions, or recommendations on how to properly structure this within the INLA/inlabru framework. I am also open to any suggestions for improving the code or model specification.

Thanks in advance

Regards

Hector Ramírez

ramire...@gmail.com

[Finn Lindgren]

Hi,

you left out the part of your code that actually defines the point process observation model, but if you've defined that correctly that shouldn't need to change.
Since you store the covariates with one layer per time index, you just need to supply the time index to the eval_spatial function.
There's no need for the "f...()" approach, as you're simply feeding the data to eval_spatial, so you can simplify the component definitions to

(assuming that "t" is a discrete time index starting at 1)

formula <- geometry + t ~ Intercept(1) +

Aspect(Aspect, main_layer = t, model = "linear") +
Curvate(`Cross-Sectional Curvature`, main_layer = t, model = "linear") +
DEM(DEM, main_layer = t, model = "linear") +
Eatness(Eatness, main_layer = t, model = "linear") +
Slope(Slope, main_layer = t, model = "linear") +
spatio_temporal_field(list(space = geometry, time = t), model = stmodel)

(Something went wrong when you copied the component definition, so I filled in the last bit...)

When you provide a terra raster as input and a specified layer, the default is to call eval_spatial(the_raster, .data., layer = main_layer),

after evaluating `main_layer` in the data context, i.e. usually it's the same as .data.$layer

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/131f1b6b-1a31-43bf-8a79-0a03cf0af2d8n%40googlegroups.com.

--

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

[Finn Lindgren]

There is bit more info in the documentation, see ?bru_comp

Finn

[H. Andrés Ramírez Maguiña]

Hello Finn,

It is a pleasure to meet you through this platform.

Thank you very much for your quick response. I just realized that I omitted the information you noticed.

Here is my complete original formula:

formula <- geometry + t ~ Intercept(1) +
Aspect(f.Aspect(.data.), model = "linear") +
Curvate(`f.Cross-Sectional Curvature`(.data.), model = "linear") +
DEM(f.DEM(.data.), model = "linear") +
Eatness(f.Eatness(.data.), model = "linear") +
Slope(f.Slope(.data.), model = "linear")

+
field( list(space = geometry, time = t), model = stmodel) )

where stmodel is defined as:

stmodel <- stModel.define(

smesh, tmesh, model,
control.priors = list(
prs = c(5, 0.05),
prt = c(3, 0.5),
psigma = c(2, 0.05)),

constr = TRUE)

and here is the lgcp() call:

fit <-lgcp(formula,
             data = dis_subset,
             domain = list(
               geometry = smesh,
               t = tmesh),
             sampler = crop_area,
             options = bru_options_set(
               verbose = TRUE,
               bru_verbose = TRUE,
               control.compute = list(
                 dic = TRUE,
                 waic = TRUE,
                 cpo = TRUE)))

where dis_subset is my dataset, a data.frame containing the coordinates (x, y), the column t (four-month periods from, let's say, 1 to 10), and the point geometries as an sf object. crop_area corresponds to the crop area where the sampling was conducted.

I followed your suggestion to simplify the component definition for the covariates that do not vary over time, and everything worked well in that case. However, I encountered some issues when incorporating the time-varying covariates.

Here is the procedure I followed. For simplicity, I tested the model using only one temporally varying covariate, for example NDVI. This is stored in a list called "rasters", and "rasters$NDVI" is a terra SpatRaster object with 10 layers. It is worth mentioning that the values in the t column of dis_subset, corresponding to each observation point, range from 1 to 10, and likewise each raster layer is named from 1 to 10.

This is the formula including the NDVI covariate:

formula <- geometry + t ~ Intercept(1) +

NDVI(rasters$NDVI, main_layer = t , model = "linear") +
field( list(space = geometry, time = t), model = stmodel) )

and the rest of the model remains exactly the same as before:

fit <- lgcp(...)

As a result, I obtained the following error message:

Error en .subset2(x, i, exact = exact): subscript out of bounds

I tried some modifications without success, which I show below together with their corresponding console messages:

formula <- geometry + t ~ Intercept(1) +

NDVI(rasters$NDVI, model = "linear") + #without main_layer = t
field( list(space = geometry, time = t), model = stmodel) )

or

formula <- geometry + t ~ Intercept(1) +

NDVI(rasters$NDVI[[t]], model = "linear") + #without main_layer = t
field( list(space = geometry, time = t), model = stmodel) )

The model starts running normally, but eventually reaches a point where unusual maxld values appear (usually these values remained between approximately -8300 and -8100 in my previous models), and the following warning message is repeated indefinitely:

maxld= 6265.1026 fn=110 theta= 3.9058 1.7003 5.4916 [1502.72, 23.414]
maxld= 6438.9219 fn=111 theta= 3.9054 1.7022 5.4984 [1518.61, 23.258]
*** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=568
*** WARNING *** GMRFLib_2order_approx: reset counter for 2 NAN/INF values in logl

Normally, my models with time-independent covariates take no more than about 15 minutes to run, but this warning continues to persist even after 40 minutes.

Finally, have tried main_selector = t, with this error:

Error in `where[[selector]]`: ! Can't extract column with `selector`. ✖ `selector` must be numeric or character, not a <standardGeneric> object. Run `rlang::last_trace()` to see where the error occurred.

I would greatly appreciate your support, and thank you in advance.

[Finn Lindgren]

I think the issue is that you’re giving tmesh as the integration domain, when it should be 1:10

Your latent model is built on continuous time, but your data uses discrete time.

It’s ok to use tmesh for the field component, but the integration domain needs to be discrete, as continuous domain integration will evaluate at points between the integers as well.

Check by calling eval_spatial() directly to verify what comes out for different values of the.

The  [[t]]     approach can’t work.

Finn

On 28 May 2026, at 23:57, H. Andrés Ramírez Maguiña <ramire...@gmail.com> wrote:



To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/27bebd9c-c477-4952-bc25-04c125268300n%40googlegroups.com.

[Finn Lindgren]

To use the main_selector argument you would need main_selector=“t”.

To check covariate evaluation,

eval_spatial(theraster, where=fm_int(domain, samplers), selector=“t”)

Where domain and samplers are the same as when calling lgcp().

You can then check if you get NAs or other problems when the t domain is tmesh, vs when it is 1:10.

Finn

On 29 May 2026, at 06:50, Finn Lindgren <finn.l...@gmail.com> wrote:

I think the issue is that you’re giving tmesh as the integration domain, when it should be 1:10