Creative methods for SPDE models in mountainous areas

[Chris Smith]

Dear R-INLA (and especially Finn who has helped me in the past),

      First, thanks again for creating these tools for dealing with spatial autocorrelation...my research would not be possible right now with out these.

I am running an INLA SPDE model on ~1000 pts scattered across roughly 150km of the northern Sierra. In a highly mountainous landscape, obviously Using a 2-D mesh has conceptual problems arise when points are close in space, but far apart in elevation (e.g. points at high elevation should have high occupancy, but they may be close spatially to pts low in elevation with low occupancy).  I have tried several things that were recommended on this chat, and wanted to document them here for other user, and get help, as my models still don't seem to be dealing with this effectively.

I have been modeling the slope of the terrain on the range parameter (steeper slopes have less autocorrelation distance/range), using this code kindly provided by Finn.  For this, I have been using a slope raster with a moving window at 250m to avoid micro-topography problems.

 

sigma0 <- 1
size <- fm_diameter(mesh)      #max diameter of mesh is 215 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))

Several lines of evidence suggest this isn't working for me. First, Theta4 (which I believe corresponds to the slope of the covariate “slope” on the autocorrelation range) is positive. From what I understand, this means that as slope of terrain increases, so does the range of the autocorrelation which is the reverse of what I want. Secondly, Theta4 is tiny and barely significant (Theta = .008, CI=-.02, .038) in my models (slope of terrain is in degrees and all rasters are projected into kilometers). Finally, when I run a stationary model with no slope covariate on the range, and WAIC are only 7 apart.

I have tried several other things, including decreasing the mesh size (from max edge of 1km down to .3 km) to better approximate the landscape, using a raster for slope at 35m pixels, and getting rid of a fixed effect in the model that is strongly correlated with elevation to see if it would change Theta4, but there is little change. I also tried incorporating a z-coordinate into my Field, but it also had little impact:

coords_3d_data <- as.matrix(data.frame(
  easting  = data$easting, northing = data$northing,  z_km     = data$talus_elev_mean / 1000 ))

 

My questions are

1. I am avoiding including a fixed effect of elevation because it is highly correlated (r=.85) with other covariates like Mean Annual Temperature that are important to investigate. Is there any solution you can think of using fixed effects?

2. Is there any clever way to use a barrier model to help make areas with very different elevations (steeper terrain) less autocorrelated? (or anything else you can think for me to try).

 3. Is there a way for me to test if my model is behaving reasonably without the slope covariate on autocorrelation? Preliminary predictions for example look reasonable.

Thanks for any help you can give.

Kind regards,

Chris