Help combining different INLA packages to model mountainous terrain

[Chris Smith]

  Dear R-INLA community

       First, thanks again for providing awesome tools for us ecologists dealing with spatial autocorrelation. I am a PhD student trying to run an SPDE model in highly mountainous terrain with a large lake (Lake Tahoe) in the middle. I am trying to get help on the combination of packages and model types to use. As I understand it, I have several options:

inlabru: makes model specification easier (but inlabru's spde.posterior() does not work for non-stationary models)

rSPDE: estimates the smoothness parameter

barrier model: models barriers more accurately than simply inserting a polygon "hole"

non-stationary models: Finn graciously provided code for modeling "slope" of the terrain              on the mesh (using kappa and tau), which should help dissasociate points close                 in space but far in elevation 2000m vs 3000m (which can be within 1km of each 

             other).

The objectives are to look at fixed effects on the occupancy of a small mammal (pika), and to plot the decline of the autocorrelation over distance/range (n=1000 data pts, study area = 150x60km).

Question 1: 

I am worried that using an SPDE barrier model (https://eliaskrainski.github.io/INLAspacetime/articles/web/barrierExample.html) will make specifying a covariates on kappa and tau quite difficult, as both are non-stationary models. I was thinking for this reason, to simply put a 'hole' in my study area for Lake Tahoe (the older method), as the correlation in steep terrain is of much higher concern. Does this seem reasonable? (I am still new to INLA;  https://becarioprecario.bitbucket.io/spde-gitbook/ch-intro.html#sec:mesh 

Question 2:

I was thinking of using rSPDE to estimate the smoothness parameter, while modeling slope as a covariate on kappa and tau. It seemed like estimating the smoothness parameter would be important in steep mountainous terrain. However, I ran a basic rSPDE model on my data, and nu = 0.7, which I believe means alpha= 1.7 or close to 2 or the default in INLA...so maybe rSPDE is not necessary?  I was unsure if rspde.result() could also be called to plot the autocorrelation over distance (range).  Any suggestions on this workflow?

https://davidbolin.github.io/rSPDE/articles/rspde_inla.html#predictions 

I am still learning all of the different packages associated with INLA and their functions, and any help understanding how to best combine these for my mountainous system would be appreciated.

Thanks for any help you can provide,

Chris

[Finn Lindgren]

Hi Chris,

For question 1, the barrier model is really just a specific special case of a non-stationary model.

The main difference between the implementation of inla.spde2.matern with non-stationary tau/kappa parameters and in barrier model implementations is that inla.spde2.matern uses per-mesh-vertex parameter values, and the barrier model implementations use per-triangle values, which makes "holes" work a bit better in the cases when the mesh edges are aligned with the edges of holes. So you can get a similar model with inla.spde2.matern by having a basis function that is the indicator for whether you're "inside the barrier", and wither let it estimate the effect of the "hole", or include it as a strong feature of the mean structure of the parameters (i.e. in the first column of B.tau and B.kappa). If the mesh has sufficient resolution, this should be very different from a barrier model.

For question 2 and non-stationary models, I'd just use alpha=2 and not estimate the smoothness.

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/b24de723-4042-463c-a273-29910ecd9285n%40googlegroups.com.

--

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

[Chris Smith]

Hey Finn,

      One more question if you don't mind on this. I have been using "slope" of the terrain as a covariate on kappa to help deal with the autocorrelation being less when elevation (slope) is changing rapidly. I have been running into some issues recently and was wondering if you have any creative thoughts on how a barrier model could be implemented to model a similar effect. The idea is to help points that are close in 2 dimensional space, but are often 1000m apart in elevation be less autocorrelated than in flat terrain.

Thanks for any help you can give,

Chris

[Chris Smith]

And apologies....one more question. Previously, I was using a non-stationary model with code you kindly helped supply:                

sigma0 <- 1

size <- fm_diameter(mesh) #max diameter of mesh is 475 km

range0 <- size/5
kappa0 <- sqrt(8)/range0
tau0 <- 1/(sqrt(4 * pi) * kappa0 * sigma0)

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

spde.slope <- inla.spde2.matern(mesh,
                B.tau = cbind(log(tau0), -1, -z1, +1, +z1),
                B.kappa = cbind(log(kappa0), 0, 0, -1, -z1),
                theta.prior.mean = rep(0, 4), theta.prior.prec = rep(1, 4))
 

In the output, I am looking at 3 variables (using Random Walk 2, similar to a GAM) to look at their effect on this animals occupancy.  The output right now I believe is showing that my variable on kappa (slope) is not having a significant effect (Theta 2 and 4). Am I interpreting this correctly?  If so, do you have any suggestions, how I might change how kappa is being specified?   

mean sd 0.025quant 0.5quant Precision for inla.group(MAT, n = 15) 8.53e+10 8.21e+13 0.000 27.918 Precision for inla.group(Hiemel, n = 15) 1.83e+08 6.58e+11 0.000 24.434 Precision for inla.group(Talus4km, n = 15) 5.65e+13 1.01e+16 0.004 442.225 Theta1 for field 1.22e+00 3.13e+00 -4.811 1.179 Theta2 for field -4.00e-03 1.76e-01 -0.359 -0.001 Theta3 for field -4.07e-01 2.98e+00 -6.147 -0.451 Theta4 for field -1.10e-02 1.60e-01 -0.333 -0.008

[Chris Smith]

[Finn Lindgren]

SOmething seems strange in those estimates; the precision parameters for MAT/Heimel/Talus4km appear mostly unconstrained; intervals virtually 0 to near infinity, with mean larger than the 97.5% quantiles.
In light of that, I wouldn't attempt to interpret the field parameters at all.

Perhaps try running the model without the field component, or run other simplified versions of the model, to figure out what's going on with those other variables, before adding back the field component.

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/5ba0fcd6-3ae8-468d-8e6b-3a34cf413bc0n%40googlegroups.com.