Can INLA be used to fit non-linear fixed effect in LGCP?

[David Dayi Li]

Hi all,

I am fitting a LGCP using inlabru with some covariate information. However, the actual covariate relationship with the rate function is not log-linear. To be more concise, assume my covariate is X, my current model has the following form:

eta = g(X; alpha) + U,

Here g is a non-linear function of X and alpha is the parameter vector determining the structure of g. U is the spatial random effect.

I am wondering if INLA can be used to fit such a model which would also provide me with inference on the non-linear parameter alpha. If so, how can I implement this in INLA?

Thank you very much!

David

[Finn Lindgren]

Hi David,

there are two types of linearity here; whether the relationship between covariate and eta is nonlinear in the covariate, and whether the eta is nonlinear in the latent Gaussian variables. It’s not entirely clear which kind of nonlinearity you’re referring to.

Nonlinearity in the latent variables can be handled in inlabru (and is one of the reasons inlabru exists) Nonlinearity in the covariate effect in the sense of e.g. a rw2 model where the covariate is the domain can be handled in both. For constructing the lgcp likelihood, inlabru is usually the easiest alternative.

If the nonlinearity is with respect to the latent variables, you can implement it in inlabru using an R expression that calls the needed nonlinear function that computes the contribution to eta. The tricky part here is the it needs to be constructed so inlabru itself can evaluate it at new locations. The details of this depends on your precise model definition.

 For things like an rw2 model of a spatial covariate, see the online covariate examples for inlabru for some syntax options.

Depending on the details on your model structure, there may be specific tricks, so I can’t say much more without a more complete mathematical model definition for the desired predictor.

Finn

On 1 Aug 2021, at 19:23, David Dayi Li <davi...@gmail.com> wrote:

Hi all,

--
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/47a6de9d-9af5-46cb-8437-94c2e469bd8an%40googlegroups.com.

[David Dayi Li]

Hi Finn,

I think the non-linearity I was refferring to is in the covariate effect. To be more precise, normally our model would be something like

eta = X*beta + U,

but my model here is something like (for illustration purpose):

eta = (X^alpha)*beta + U.

So the konwn covariate is related to eta through a non-linear function and I also wish to conduct inference on alpha here as well. Hope that clears things up.

Thanks,

David

[Finn Lindgren]

Ok!

For regular georeferenced observations, that would be as simple as defining alpha as some suitable transformation alpha_fun() of a N(0,1) variable, and writing the predictor formula as X^alpha_fun(alpha_component)

For lgcp models, it’s a bit harder. The alpha part would work exactly the same, but the X-evaluation becomes the harder part. (X*beta)^alpha would be easier, if you already know the sign of the effect, since then you can define X*beta as a regular linear effect of the covariate, and inlabru knows how to spatially evaluate those. But for the model as you wrote it, you have (at least) two options:

1. Define a function that evaluates the spatial covariate at the input locations and call that function in the predictor formula, or

2. Define a linear predictor component for X, but set the latent mean to 1 and the precision to a large value.

The kind of function needed for X in alternative 1. Is the same type as function used in some of the covariate examples for inlabru (there is an bru_fill_missing() function that fills missing values with the nearest available value, which is usually preferable to the fill_with_zero used in the old examples.)

Finn

On 1 Aug 2021, at 19:49, David Dayi Li <davi...@gmail.com> wrote:

Hi Finn,

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/ee9ab1b2-2c28-4d98-bcae-869c922c7276n%40googlegroups.com.

[David Dayi Li]

Hi Finn,

Since there does not seem to be any code example available for this, I am not exactly sure how to implement this. Also, I am not sure I understood some of the points you mentioned.

First off, is the N(0,1) variable for alpha_fun() a prior distribution for alpha? If it is, then I don't think this is suitable since I know alpha is strictly positive. Also, I do know that the sign of beta is positive, but I am not sure what you mean by defining X*beta as a linear effect of the covariate. 

Would you be able to provide or point me to a code example of how to imlement this? Thank you very much!

David

[Finn Lindgren]

Hi David,

since the predictor expression in inlabru in principle can be any R expression, it's difficult to cover all possible models in the examples, especially those that are essentially current research problems (in terms of how well the iterative method and the resulting posterior approximation behaves, etc).

But with your additional details, here's how one can do it in your case:

alpha_fun can be used to specify a specific prior distribution, yes. Since you say alpha is known to be strictly positive, that just means you need a strictly possible transformation there.

For example, say you want an exponential distribution for alpha. Then qexp(pnorm(u),rate=...) has an exponential distribution with given rate parameter, if u is N(0,1). Replace qexp() with the quantile function for whatever prior distribution you want for alpha.

Below I'll assume X is a SpatialGridDataFrame where the first data column is the covariate.

alpha_fun <- function(u) {

  qexp(pnorm(u), rate=prior_rate)

}

comp <- ~ - 1 + Xbeta(X, model = "linear") +

   alpha_internal(1, model="linear", mean.linear=0, prec.linear=1) +

   U(coordinates, model=spde)

form <- coordinates ~ abs(Xbeta)^alpha_fun(alpha_internal) + U

The model above defines eta = |X*beta|^alpha + U. The absolute value is just there to guard against X*beta becoming negative, which would break when alpha becomes a non-integer.

In order for the iterative method to work, you need to initialise the beta coefficient (so the derivatives aren't zero in the first iteration):

fit <- bru(components = comp, like(model="cp", data=..., samplers=...), options=list(bru_initial=list(Xbeta=1)))

Note:  The alpha_fun definition can be made more numerically robust by using branching for pos/neg u, and log.p and lower.tail arguments to qexp and pnorm. I plan to add automated helpers for this to inlabru soon, since it's such a common thing to want to do.

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/49f3f720-d100-4f07-b4d4-0e21223efae5n%40googlegroups.com.

--

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

[Finn Lindgren]

In principle I think one can also do the X^alpha*beta + U version of the model, by using that model="offset" in inlabru can be used to evaluate constants:

comp <- ~ - 1 + x(X, model = "offset") + beta(1, model="linear") +

   alpha_internal(1, model="linear", mean.linear=0, prec.linear=1) +

   U(coordinates, model=spde)

form <- coordinates ~ beta * x^alpha_fun(alpha_internal) + U

This "offset" feature should really be called "const", and may be renamed to that in the future; the name "offset" is just because it was initially just the internal implementation of offset() in inlabru.

I haven't used this feature myself yet, so please let me know if you try it, and if it works or not! (it is tested by the package tests, but those can't cover all use cases).

(Also note I didn't include an intercept parameter in the code examples; add that in if you need one...)

Finn

[David Dayi Li]

Hi Finn,

Thanks so much for this. However, I just tried your code, I got the following error:

Error in like(model = "cp", data = spv, samplers = region, domain = list(coordinates = v_mesh)) : 

  unused argument (model = "cp")

What went wrong here?

David

[David Dayi Li]

Hi Finn,

It works now. Thanks so much for your help!

David

On Sunday, August 1, 2021 at 4:58:10 PM UTC-4 finn.l...@gmail.com wrote:

[David Dayi Li]

Also, a final question. The inferred alpha_internal parameter should not be the actual alpha parameter right? The actual distribution of alpha should instead be the one that is obtained after transforming alpha_internal using alpha_fun(), is that correct?

David

On Sunday, August 1, 2021 at 4:58:10 PM UTC-4 finn.l...@gmail.com wrote:

[Finn Lindgren]

Great!

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/9ea8f837-f3a3-4c37-ad36-95202aadf923n%40googlegroups.com.

[Finn Lindgren]

Yes, that's correct.

The alpha_internal-->alpha transformation can be seen as a 1-dimensional copula transformation; In the prior, alpha_internal is N(0,1) and alpha has the distribution determined by the chosen quantile function. In the posterior, alpha_internal has some posterior distribution, and alpha has the distribution obtained when feeding that distribution through the quantile(pnorm()) transformation.  The inla.tmarginal() function from INLA can be used to calculate the posterior density for alpha. Alternatively, use the predict() method to estimate the posterior summary statistics using posterior sampling.

Here's the theory for it:

P(alpha <= a) = P(alpha_fun(u) <= a) = P(quantile_fun(pnorm(u)) <= a)

  = P(pnorm(u) <= cdf_fun(a)) = P(u <= qnorm(cdf_fun(a)))

where cdf_fun is the CDF function for the prior model for alpha and u is alpha_internal.

INLA has functions for computing the transformation from alpha_internal: Using inla.tmarginal(alpha_fun, marginal) should produce the posterior density for alpha itself, when given the marginal density output for alpha_internal.

Or inla.emarginal(alpha_fun, marginal) for just the posterior expectation.

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/739b6900-3455-40f8-bf2d-d16ed6ee0b7dn%40googlegroups.com.

[Finn Lindgren]

...and

  alpha_fun(inla.qmarginal(c(0.025, 0.5, 0.975), marginal))

for posterior quantiles.

Finn