Graphing autocorrelation covariate effects of kappa/range (INLA SPDE)

[Chris Smith]

Dear INLA community,

      I am a PhD student trying to use INLA to model the presence/absence of a small mammal, and at what distance habitat patches stop influencing one another. Simply looking at my data, it is clear that even out to 5km, I am seeing high correlation in my data and I have started using an INLA SPDE model to try to answer this question. 

Me and my advisor are trying to use the autocorrelation parameters (kappa and tau) as a metric of how the effect of an autocorrelation covariate fades to 0 over distance--basically using the spatial autocorrelation as a metric of how the connectivity of patches drops off with distance. I have been looking online, and have failed thus far to find a straightforward worked example of how someone puts a covariate on kappa and tau, and then graphs how autocorrelation changes. We are really interested in how the range of the autocorrelation (associated with kappa) changes with distance from a point, but it seems like the relationship between range and kappa is complicated.

I've tried reading this chapter, but haven't been able to follow all that is going on:

https://becarioprecario.bitbucket.io/spde-gitbook/ch-nonstationarity.html   

I was wondering if:

1. Anyone has used spatial autocorrelation as a metric of how correlation between habitat patches changes with distance (as an index of when 1 habitat patch stop influencing other patches around it)

2. Someone could share a simple example of how to model a covariate on kappa, and then create a graph from the model predictions (using theta?) of how the autocorrelation range declines with distance. 

Thanks for any help you can provide! I imagine parts of this question may be confusing, so let me know what other details I need to provide.

[Finn Lindgren]

Hi,

for non-stationary fields, to get a parameterisation that's more interpretable in terms of correlation range and variance, see around page 5 of the Lindgren and Rue 2015 JSS paper:

  https://www.jstatsoft.org/article/view/v063i19

The model construction works the same in the current version of the package, but you should use inlabru to do the estimation to take advantage of the easier interface for spatial models;

instead of "+ f(field_index, model=spde_object) you'd use "+ field(geometry, model=spde_object)" if you have data stored in sf format, for example, and point process likelihood support is built in.

See the top menu of https://inlabru-org.github.io/inlabru/ for many examples.

The paper shows how to have the covariates affect the log-range, instead of the log-kappa values (kappa and range are directly linked, but kappa also affects the variance; the paper shows how to parameterise it using the range in a way that's decoupled from the variance). (the newer pcmatern model implementation doesn't have direct non-stationary support, but uses this technique internally for stationary models)

The INLA spde.result() function takes an estimated model and computes pointwise parameter estimates, so you can see/plot how the range parameter varies across space (this may be sufficient for your purposes).

For the shape of the covariance, the only real way to do that for nonstationary fields is to compute the covariance matrix for a collection of points.

For stationary fields, inlabru has the spde.posterior() method that can help; see "?inlabru::spde.posterior()"; it doesn't work for non-stationary models, as it assumes a stationary covariance form.
For non-stationary, one needs to construct the FEM-matrices with fm_fem() and construct the precision matrix from that; there may be partial helper functions for this out there, but we don't yet have a clear interface for that for the non-stationary case; fmesher::fm_matern_precision() only handles the stationary case.

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/c9f102b2-c1ed-40fa-b96e-75ed282f20bfn%40googlegroups.com.

--

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

[Finn Lindgren]

Correction: The spde.result function in INLA is called inla.spde.result, see "?INLA::inla.spde.result"

Finn

[Chris Smith]

Thanks Finn. I really appreciate how fast you and Havard are at returning messages on this board. It makes using INLA easier for us beginners.  I found this chunk of code in the supplementary materials, and assume this is close to what you are talking about?

sigma0 <- 1
size <- min(c(diff(range(mesh$loc[, 1])), diff(range(mesh$loc[, 2]))))
range0 <- size/5
kappa0 <- sqrt(8)/range0
tau0 <- 1/(sqrt(4 * pi) * kappa0 * sigma0)
spde <- inla.spde2.matern(mesh, B.tau = cbind(log(tau0), -1, +1), B.kappa = cbind(log(kappa0),
  0, -1), theta.prior.mean = c(0, 0), theta.prior.prec = c(0.1, 1))

###################################################
truevar <- (3 * sigma0)^2
truerange <- range0


If you have time for one more question:

Is there any built in method in INLA for calculating the amount of autocorrelation left in an INLA model after it has run? I'm thinking of extracting residuals and calculating a Moran's I value, but wasn't sure if there was a more efficient method.

Thanks again,
Chris

[Finn Lindgren]

That's the code for the stationary version yes. (In the stationary case, it's better to use inla.spde2.pcmatern(), that already uses range and sigma as parameters).

If you have covariates z1 z2 evaluated at the mesh nodes, you can use

sigma0 <- 1
size <- fm_diameter(mesh)

range0 <- size/5
kappa0 <- sqrt(8)/range0
tau0 <- 1/(sqrt(4 * pi) * kappa0 * sigma0)
spde <- inla.spde2.matern(mesh,

  B.tau = cbind(log(tau0), -1, -z1, -z2, +1, +z1, +z2),
  B.kappa = cbind(log(kappa0), 0, 0, 0, -1, -z1, -z2),
  theta.prior.mean = rep(0, 6), theta.prior.prec = rep(1, 6))


If you have covariates available in a terra raster, and they are smooth enough that the pointwise values at the mesh nodes give a sensible smooth approximation,
you can get them via the inlabru function eval_spatial;

z1 <- eval_spatial(covariates_in_terra_raster, where = fm_as_sfc(mesh, format = "loc"), layer = "z1")
z2 <- eval_spatial(covariates_in_terra_raster, where = fm_as_sfc(mesh, format = "loc"), layer = "z2")

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e54ae0c1-915a-44f2-9981-f0c0d7bfd077n%40googlegroups.com.

[Finn Lindgren]

For the last question,

> Is there any built in method in INLA for calculating the amount of autocorrelation left in an INLA model after it has run? I'm thinking of extracting residuals and calculating a Moran's I value, but wasn't sure if there was a more efficient method.

Nothing built-in, as it depends so much on the application. Some kind of residual analysis would make sense, yes.
The inlabru generate() and predict() methods can evaluate general R expressions for each posterior sample, so what you can compute is mainly up to your imagination...

In the case of no link function,

predict(fit, newdata = test_data, formula = ~ the_response_variable - (the + sum + of + the + model + component effects))

would give you the individual posterior expectation and std.dev of the residuals; E(y_test - \mu | y_estimation) and sqrt of Var(y_test - \mu | y_estimation)

(Here, \mu is the predictor expression, and the samples the calculations are based on are the posterior samples of \mu)

One can also extract some per-estimation-observation information from the fitted model object, but for general models I find it easier to use the more flexible posterior sampling approach.

Finn