Develop and test multi-scale model

[James-Thorson-NOAA]

It is sometimes useful to have capability to model spatial correlations arising on radically different scales. Examples include models for transect data, where spatial correlations within vs. among transects occur at different orders of magnitude. This can be seen e.g., here

library(INLA) 
library(RandomFields)

# Sample size
n_clust = 30
n_within = 30
sigma_within = 0.001

# RF
sd_fine = 0.5
sd_coarse = 1
SpatialSimModel = "Matern"
Nu = 1
range_fine = sigma_within
range_coarse = 0.2

# Sampling
Mu = 10

# derived
kappa_fine = sqrt(8*Nu) / range_fine
kappa_coarse = sqrt(8*Nu) / range_coarse

# Simulate locations
loc_iz = matrix(runif(n_clust*2),ncol=2)
loc_iz = loc_iz[ rep(1:n_clust,each=n_within), ]
loc_iz = loc_iz + sigma_within*matrix(rnorm(n_clust*n_within*2),ncol=2)

# Simulate RF
RM_fine = RMmatern(nu=Nu, var=sd_fine^2, scale=range_fine*2)
RM_coarse = RMmatern(nu=Nu, var=sd_coarse^2, scale=range_coarse*2)
Omega_fine = RFsimulate(model=RM_fine, x=loc_iz[,1], y=loc_iz[,2])@data[,1]
Omega_coarse = RFsimulate(model=RM_coarse, x=loc_iz[,1], y=loc_iz[,2])@data[,1]
Omega = Omega_fine + Omega_coarse

# Simulate data
lambda = exp( Omega + Mu )
n_i = rpois( n=length(lambda), lambda )

# create mesh
mesh = inla.mesh.create( loc_iz, plot.delay=NULL, refine=TRUE)

It is possible to model this appropriate using an additive function of two GMRFs (similar to this generative model), and this could be done in VAST by changing log_kappa1 and log_kappa2 from a scalar to a vector that is constant by default, with one value per Omega/Epsilon/Xi/Phi term, and using the formula interface and spatially-varying responses to estimate additional GMRFs that have the action of Omega or an IID Epsilon term.

I will this this as a development topic, and would implement if it someone volunteers to help me explore it.