SPDE Spatial Autocorrelation in areas with high elevation change

[Chris Smith]

Dear INLA google group,

     First, thanks for making such an awesome tool for ecologists out there who struggle with spatial autocorrelation in our studies.  This package recently solved a huge problem I was having with spatial autocorrelation in my data.

My study organism (American Pika) live in mountainous areas with large changes in elevation. For example, within a 1 km distance, I go from 6000 to 10000 ft up the side of a mountain and pika occupancy goes from near 0% probability up to 100%.  My understanding is that SPDE models using a mesh to model autocorrelation assume points that are closer in horizontal distances are more similar. The problem for pika is that for 2 points 1 km horizontally apart (both at 6000ft elevation), this works well. However for 2 points 1km apart at 6000ft and 10000ft, the correlation is extremely different (vs 1 km horizontally). 

I have run preliminary models with the SPDE random effect, while including elevation as a fixed effect, and the predictions seem reasonable.

I wanted to ask:

1) Does adding elevation as a fixed effect help overcome this problem, with horizontal vs vertical spatial autocorrelation?

2)  Are there other tools in INLA I should be using for 3-D autocorrelation?

3) I am worried about putting highly correlated variables like snowpack into models with elevation; do you have any suggestions how to still include correlated climate variables if I am putting elevation in all my models?

I hope this post can be helpful for other researchers working in mountainous ecosystems where elevation change is a major part of our system.

Thanks for any help you can provide.

[Finn Lindgren]

Hi Chris,

1)
there are two aspects here that are often confused; coherently spatially changing mean value, and spatial correlation; with just one realisation of the field, one can't really tell those apart, and one has to use ones modelling choices to get the kind of information one is interested in. Sometimes, it's enough to include the covariate in the predictor (i.e. controlling the mean value), but if one truly needs it to control the dependency structure (i.e. to _decouple_ the values at high and low elevation) one needs to enter it as covariates in the dependence structure. In your case, you likely need both.

The inla.sdpe2.matern() model constructor allows you to include covariates evaluated at the mesh nodes, e.g. if elevation is available as a terra raster:

library(inlabru)

elev_on_mesh <- eval_spatial(elevation, mesh$loc)

spde <- inla.spde2.matern(mesh, alpha = 2, B.tau = cbind(0, 1, elev_on_mesh, 0, 0), B.kappa = cbind(0, 0, 0, 1, elev_on_mesh))

This lets both tau and kappa vary across space. Also see https://www.jstatsoft.org/article/view/v063i19 for more details on the parameterisation, and how to reparameterise the model to be more interpretable, in terms of log(std.dev.) and log(range) instead of the less directly interpretable tau and kappa parameters.

The above is just a simple example of how to write the code. In your problem, it's more likely that you want kappa to depend on the magnitude of the _slope_ of the elevation field, so that a large slope allows it to locally have a short correlation length, which has the effect of decoupling areas of low and high elevation. (this would be a "soft" version of the "barrier" model, which is discussed elsewhere)

2) There is support for full 3D tetrahedralisations (fmesher::fm_mesh_3d()), but that is unlikely to be helpful in your case, and is more suited for e.g. 3D seismology etc.

3) If the covariates are highly correlated, then you'll get confounding in the estimates, yes. There's no universal solution to that.
You will be ba able to see in the estimates if it is a problem, e.g. if the effect of each individual component has a large posterior predictive standard deviation, but the combined predictions have reasonable std.ev.; that would mean that the estimates are negatively correlation, which is what one would expect in a confounded situation. (This often happens between the intercept and spde components as well; in those cases, just focus on their combined effect instead of attempting to interpret the individual effects).

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/d51e8c46-b6bb-4aa9-9200-a13e8d626a4en%40googlegroups.com.

--

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

[Chris Smith]

Dear Finn,

     Thanks so much for this. This has very much changed our analysis and made it much more realistic, knowing that we can now model the autocorrelation off of covariates. 

I have run some preliminary models using slope, as I agree that using elevation to model the autocorrelation doesn't make much sense (e.g. a high elevation has high autocorrelation and a low one has low autocorrelation doesn't fit our system). I tried running a model using elevation (scaled) as a fixed effect (elevation in our system is a reasonable predictor of pika) and slope in degrees to model spatial autocorrelation. I then plotted the results and compared them to a null autocorrelation model with elevation (scaled) as a fixed effect and no covariates on the autocorrelation. The good news is both results look quite reasonable, but the bad news is that they look identical :).

It seems to me that slope is not having almost any effect on the predictions. I have printed out the predictions of both models below.

Model with no spatial autocorrelation covariates: 

Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld b0 -2.663 1.182 -5.140 -2.602 -0.505 -2.416 0

elev 2.070 0.315 1.474 2.062 2.715 2.063 0

mean sd 0.025quant 0.5quant 0.975quant mode

Theta1 for s 6.508117 0.2514101 6.067679 6.493699 7.050241 6.416033 Theta2 for s -10.222326 0.9315390 -12.368137 -10.102403 -8.814396 -9.501286

Model with slope autocorrelation covariate:

mean sd 0.025quant 0.5quant 0.975quant mode kld b0 -3.582 1.707 -7.281 -3.418 -0.668 -2.846 0 elev 2.029 0.298 1.445 2.028 2.616 2.028 0  

mean sd 0.025quant 0.5quant 0.975quant mode Theta1 for s 6.650 0.033 6.576 6.653 6.703 6.669 Theta2 for s -0.012 0.009 -0.031 -0.012 0.005 -0.010 Theta3 for s -9.611 0.030 -9.664 -9.613 -9.547 -9.622 Theta4 for s 0.022 0.010 0.005 0.022 0.042 0.021

Do you think this suggests I need to find a different parameter to try to model autocorrelation?

I also wondered about the scale of the mesh compared with the scale of the slope variable. We are using UTM's and currently have a max edge of 2000 with ~8000 mesh points. I noticed that B.tau and B.kappa here seem to be extracting slope at each of the mesh nodes. Given that many of our valleys may be less than 1km wide with significant elevation/slope change in the middle, should I be setting max edge to be a smaller value so the mesh picks up these smaller scale elevation/slope changes?

Thanks again for any help you can give, and let me know if you need more information. I hope this can be useful for other mountain researchers facing similar challenges.

[Chris Smith]

Dear Finn,

    My sincere apologies that this output printed so unreadibly. It was fine in my word document, but not when it posted in the email.  I have pasted it below so this hopefully is easier to read.

Kind regards and thanks again for all your help,

Chris

Model with no autocorrelation covariates:

Fixed effects:

         Mean        sd      0.025quant    0.5quant       0.975quant    mode    kld

b0    -2.663   1.182   -5.140               -2.602            -0.505           -2.416  0

elev  2.070    0.315   1.474                 2.062             2.715             2.063  0

 

                          mean       sd                 0.025quant    0.5quant    0.975quant   mode

Theta1 for s  6.508117   0.2514101   6.067679      6.493699    7.050241      6.416033

Theta2 for s  -10.222      0.9315390  -12.368137   -10.102403 -8.814396    -9.501286

Model with slope autocorrelation covariate:

           mean    sd            0.025quant    0.5quant  0.975quant    mode kld

b0     -3.582      1.707     -7.281            -3.418       -0.668             -2.846 0

elev   2.029       0.298      1.445             2.028        2.616               2.028 0  

 

                       mean      sd        0.025quant     0.5quant 0.975quant    mode

Theta1 for s  6.650     0.033     6.576               6.653       6.703            6.669

Theta2 for s  -0.012    0.009    -0.031             -0.012       0.005           -0.010

Theta3 for s  -9.611    0.030    -9.664             -9.613      -9.547           -9.622

Theta4 for s    0.022   0.010      0.005              0.022       0.042             0.021