how to define functions of covariates in lgcp

[ASeatonSpatial]

For a straightforward e.g. Poisson glm say, I can include any covariates with the data passed to bru(). I can then define functions to transform those values in my predictor:

my_fun = function(x, theta){
   # something involving x and theta
   }

# covariates stored in column z in the data passed to bru()
fml = y ~ my_fun(z, theta) 
cmp = ~ theta

How can I do something like this when the model is a log-Gaussian Cox process? In this case the covariate is defined everywhere, not just at the observed points.

I am considering something like

my_fun = function(x,y, theta){
   # get covariate value at location x,y
   # do something with this value that depends on theta
   }

This is similar to how do continuous covariates

my_cov = function(x,y)  { # get cov value at x,y}
cmp = ~ cov(main = my_cov, model = "linear")

I don't understand how the model knows what to do with this component. It assumes it needs to get coordinates (probably a 2 column matrix?), split it into an x, y arguments and call this function? I think because I don't understand this I am struggling to define some non-linear function of the covariate. How can I do something like this for a non-linear function of the spatial covariate within a lgcp model?

[finnlindgren]

For full control over component inputs, you can use the full function call (if x and y are variable names in the input data, or the coordinate names):
cmp = ~ cov(my_cov(x, y), model = "linear")
The function must be vectorised, so that it returns the results for each x[i], y[i] pair.

This will cause inlabru to first evaluate the function call and store the result as a new internal variable. This is then used in the same way as if it had been part of the input data, so for regular models this doesn't really provide anything extra; you could just as well have created that variable before calling inlabru. But for LGCP models, this approach allows it to call your evaluator function based on new data and new locations generated internally.

To use such a spatial covariate with a random effect:
cmp = ~ cov(my_cov(x, y), model = "rw2")
The model formula can then use the defined components with a general (vectorised) R expression:
form = coordinates ~ log(1 + cov^2)
The default formula is defined as the sum of all the defined model components, but it's better practice to explicitly define the desired model expression.

The shorthand main = my_fun should perhaps be avoided, since it's not explicit about what will be given to the function; it should be equivalent to main = my_fun(.data.), which gives it the whole data object. For SpatialPointsDataFrame data input, you can choose to have your function take a coordinate matrix as input, and use main = my_fun(coordinates(.data.)). This is also effectively how the shorthand main = coordinates works; it's a function (from the sp package) that returns the coordinate matrix when called with a SpatialPoints. However, the special case "coordinates" is also detected and treated in special ways.

[ASeatonSpatial]

Thanks Finn this is clearer to me I think. Let's see if I understand. For the lgcp this .data. object is actually the combined integration scheme locations and observed point locations, the Berman-Turner type thing that is in the "Going off grid" paper. This is a spatial points dataframe (at this point within the internals of inlabru anyway). So an alternative is to not do a vectorised form and to just work with this whole data object.

Then one uses this function to define a model component with model = "linear". In the model formula we can do other things like functions of this component.

One question that comes to mind is usually a model component with model = "linear" implies a coefficient associated with this component. Is that correct? I might not want to have that. E.g. an application I have in mind is modelling a covariate as a detection probability so I want a function that maps to [0,1]. Something like

cmp =  ~ cov(main = my_cov, map = "linear") + mu + lsig
fml = ~ pnorm((cov - mu)/exp(lsig), log.p = TRUE)

In this case I don't want any other parameters other than mu and lsig.

[finnlindgren]

Yes, .data. is the combined internal data object for the likelihood. If a model has more than one like() object, the function will be called once for each of the likelihood objects.

If the function takes the whole data object as input is still has to be vectorised to return the covariate value for each row of the data
frame.
After the internal data has been constructed and all the component inputs have been evaluated, each component works in the same way as in regular INLA, but with the added option of specifying how the inputs should be mapped to the latent Gaussian variables, via a mapper specification (and optionally also group_mapper and replicate_mapper).

In the model formula, one can use both the effect of the component by using it's label (cov) or directly the latent variables themselves (cov_latent, if like(..., allow_latent=TRUE) is used). The latter will be needed for more advanced model constructions, and can be particularly useful in predict() expressions, but usually it's the effects that's needed. In your example, you should use mu(1)+lsig(1) in the component definition, to ensure it doesn't accidentally turn mu and lsig into linear effect covariates based on some variables with those names that R might expose to the expression evaluator. +mu(1) leads to mu containing the effect vector rep(1, nrow(.data.)) * mu_latent, where mu_latent is the scalar latent variable. If the expression is defined in suitable vectorised form, either mu or mu_latent can be used, so after changing mu+lsig to mu(1)+lsig(1) in your component definition, your example should work, since pnorm is vectorised in the required way.

(Note: the mu+lsig syntax is meant to give a warning in inlabru version 2.2.5, but I'm not sure that warning always appears.)

[finnlindgren]

Just noted an error in the example; map="linear" should be model="linear". Also, main is now guaranteed to be the first parameter, so it doesn't need to be named.

[finnlindgren]

If I can figure out how some of the rlang package features work, I will deprecate .data. in favour of .data to make the expression use syntax more similar to the tidy evaluation syntax. It's not actually using any of those features, but it's a package checking issue.

[rita86ribeiro]

Hi everyone,

I would like to define a nonlinear detection function based on the nonlinear effect of a covariate in a lGCP. This covariate is a spatial pixel data frame and is the time travelling to some health offices (i.e., each grid cell provides an estimation of the time needed to travel to the health offices). Here, the detectability of the point process will decay as a sigmoid shape in relation to the time travelling. I have started to define it a half normal function, however inlabru run results in an error: Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :arguments imply differing number of rows: 371, 209330

I have:

effort <- access_senasa$effort ##access_senasa is my spatial pixel data frame and effort is the column with the data

half_normal <- function(effort, lsig) { exp(-0.5*(effort/exp(lsig))^2)}

cmp1 <- ~ Intercept(1) + sig(1)
fml1 <- coordinates ~ Intercept + half_normal(effort, lsig)

fit1 <- lgcp(components = cmp1, data=vbr, samplers=expval, domain=mesh1, formula=fml1, options=list(verbose=TRUE)

##vbr is a spatial points DF, samplers is a spatial polygon DF

In the end, I would like to reproduce this simulated example but in the LGCP framework:
n=100
beta0=5 ## mean abundance
beta1=1.5
covar=rnorm(n) ##this is another covariate part of the linear predictor
dist=runif(n,0,10)

half_normal <- function(dist,covar, beta1) {
linpred= beta1*covar
exp(-0.5 * (dist / exp(linpred))^2)}

simD1 = data.frame(dist=dist,
y= beta0 + half_normal(dist=dist,covar=covar,beta1 = beta1) + rnorm(n,0,.01))
cmp1 = ~ beta1(1) + Intercept(1)
fml1 = y ~ Intercept+ half_normal(dist,covar,beta1)
m1 = bru(components = cmp1, formula = fml1, data = simD1, family ="gaussian", options = list(verbose=T))

Thank you

sessionInfo:
R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS 13.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base

other attached packages:
[1] rgeos_0.5-9 readxl_1.4.0 dplyr_1.0.8 ggpolypath_0.1.0 ggpubr_0.4.0 RColorBrewer_1.1-3 rgdal_1.5-32
[8] ggplot2_3.3.6 INLA_22.04.09 foreach_1.5.2 Matrix_1.5-1 inlabru_2.6.0.9003 raster_3.5-15 sp_1.5-0

[finnlindgren]

Hi, since for lgcp model, inlabru needs to evaluate the covariate at various locations, the code needs to access those location in the formula. for direct component inputs, gridded covariates are auto-evaluated, but for uses like yours you have to explicitly call the evaluation function (also, since the predictor is on log-scale, you need to add the log of the detection probability, e.g. as below):

effort <- access_senasa$effort ##access_senasa is my spatial pixel data frame and effort is the column with the data
log_half_normal <- function(effort, lsig) { (-0.5*(effort/exp(lsig))^2)}
cmp1 <- ~ Intercept(1) + sig(1)
fml1 <- coordinates ~ Intercept + log_half_normal(eval_spatial(access_senasa, .data., layer = "effort"), lsig)
fit1 <- lgcp(components = cmp1, data=vbr, samplers=expval, domain=mesh1, formula=fml1, options=list(verbose=TRUE)

For a fully working eval_spatial method, you may need a more recent inlabru version; see https://inlabru-org.github.io/inlabru/news/index.html for info about when some features were introduced and/or had bug fixes.

[rita86ribeiro]

Thank you so much Finn for your help. This is great.

[osupplisson]

Hi,

I have a related question and have not been able to figure out the right syntax.

I currently have a model running fine, essentially specified as :
$y_j(s,a) \vert (p_j(s,a),N(s,a)) \sim B(p_j(s,a),N(s,a), j \in \lbrace 1, 2\rbrace$
$p_j(s,a) = plogis(b_0 + f(a) + f_j(a)+W(s) + W_j(s))$

with j a phenotype, a the age, and s the 2D spatial domain. All W and f are specified through SPDE.

Following joint models for preferential sampling (e.g. #57), I'd like to add an LGCP component with this form:

$log(\lambda(a,s)) = c_0 + \psi(age) + \beta * W(s)$, with $\beta$ the scaling parameter for W(s) and $\psi(.)$ an SPDE

I am stuck at specifying $\psi(age)$.

I got from this post and another one from the R-INLA discussion group that I must provide a component like:

psi(main = age(x, y), map = spdeage)

with age(x, y) a function providing to inlabru the values for age at any point in space for each data point (so a function just returning the value for age from the dataset whatever the coordinates).

I haven't been able to find the right syntax using sf object (https://inlabru-org.github.io/inlabru/articles/components_vignette.html only gives example for sp objects).

Any help would be appreciated !

Best,
Olivier

[finnlindgren]

Hi, the components_vignette document should be long gone from the webpage! It seems like pkgdown hasn't cleared away old versions, and I had missed that one in the list of redirects to renamed documents. The current version is https://inlabru-org.github.io/inlabru/articles/component.html which does mention sf objects. The short story for sf is that you change coordinates to geometry, and most of the rest remains the same. Most of the vignettes have been updated to use sf now, so please go to https://inlabru-org.github.io/inlabru/index.html and follow the links in the menu from there, so you're guaranteed to see the latest versions (I'll also add a redirect for the old one you found).

The syntax you need depends on what the model is meant to be. I don't think having age be a function of space is really what you want here, but rather that age is a property of the individuals (i.e. not a purely spatial quantity). For that, you need to define the product domain of space and age:

mesh_space <- ...
mesh_age <- fm_mesh_1d(knotsforagemodel)
cmp <- ~ 0 + ... + psi(age, model=spdeage) + W(geometry, model = spdeW) + betaW(1)
formula_lambda <- geometry + age ~ ... + psi + beta_W * W
fit <- bru(cmp,
  like(...),
  like(formula_lambda,
    family="cp",
    data = yourdata,
    domain = list(geometry=mesh_space, age=mesh_age),
    samplers = spatial_domain_polygon),
    options = list(bru_initial = list(beta_W = 1, W = 1)
)

Note: If this is not the type of model you're intending, then lambda should strictly be a function of s only (and any involvement of age would be in terms of age(s), which then would need a function to generate the values). I'd need more precise information to say more in that case.

Note 2: In most INLA examples with things like beta_W * W, you'll likely see the "copy" feature being used, to incorporate beta_w as a "hyperparameter" instead of a latent Gaussian; I've recently noticed a practical numerical problem with that, and the above alternative non-linear predictor approach currently seems to work better. (But you can try the "copy" option; in my case, it was pretty obvious when something was wrong, so if it seems to generate sensible results you're probably ok using it)

[osupplisson]

Thank you very much for your answer and example, which is indeed exactly my use case. Thanks also for the clarification about function of space vs property of the unit.

The Note is also incredibly useful because I believe I encountered a computational issue yesterday, following the inclusion of a scaled copy in a model that already had a fairly large number of hyperparameters.

Cheers,
Olivier