This vignette details the different covariance structures available in clustTMB.
Covariance | Notation | No..of.Parameters | Data.requirements |
---|---|---|---|
Spatial GMRF | gmrf | 2 | spatial coordinates |
AR(1) | ar1 | 2 | unit spaced levels |
Rank Reduction | rr(random = H) | JH - (H(H-1))/2 | |
Spatial Rank Reduction | rr(spatial = H) | 1 + JH - (H(H-1))/2 | spatial coordinates |
clustTMB fits spatial random effects using a Gaussian Markov Random Field (GMRF). The precision matrix, , of the GMRF is the inverse of a Matern covariance function and takes two parameters: 1) , which is the spatial decay parameter and a scaled function of the spatial range, , the distance at which two locations are considered independent; and 2) , which is a function of and the marginal spatial variance :
The precision matrix is approximated following the SPDE-FEM approach [@Lindgren2011], where a constrained Delaunay triangulation network is used to discretize the spatial extent in order to determine a GMRF for a set of irregularly spaced locations, i$.
Prior to fitting a spatial cluster model with clustTMB, users need to set up the constrained Delaunay Triangulation network using the R package, fmesher. This package provides a CRAN distributed collection of mesh functions developed for the package, R-INLA. For guidance on setting up an appropriate mesh, see Triangulation details and examples and Tools for mesh assessment from
In this example, the following mesh specifications were used:
loc <- meuse[, 1:2]
Bnd <- fmesher::fm_nonconvex_hull(as.matrix(loc), convex = 200)
meuse.mesh <- fmesher::fm_mesh_2d(as.matrix(loc),
max.edge = c(300, 1000),
boundary = Bnd
)
## Loading required namespace: INLA
Coordinates are converted to a spatial point dataframe and read into the clustTMB model, along with the mesh, using the spatial.list argument. The gating formula is specified using the gmrf() command:
Loc <- sf::st_as_sf(loc, coords = c("x", "y"))
mod <- clustTMB(
response = meuse[, 3:6],
family = lognormal(link = "identity"),
gatingformula = ~ gmrf(0 + 1 | loc),
G = 4, covariance.structure = "VVV",
spatial.list = list(loc = Loc, mesh = meuse.mesh)
)
## intercept removed from gatingformula
## when random effects specified
## spatial projection is turned off. Need to provide locations in projection.list$grid.df for spatial predictions
Models are optimized with nlminb(), model results can be viewed with nlminb commands:
# Estimated fixed parameters
mod$opt$par
## betag betag betag betad betad betad betad
## 0.1778517 0.5710441 0.1653750 2.0157772 4.3160911 5.4259797 6.7095829
## betad betad betad betad betad betad betad
## 1.0064936 3.6030483 5.2113125 6.2155353 0.1259828 3.1475076 4.2016919
## betad betad betad betad betad theta theta
## 5.2523416 -1.4361461 3.1133028 4.2118636 5.1996631 -1.2100901 -2.9055603
## theta theta theta theta theta theta theta
## -1.2794783 -1.2502294 -2.5624050 -3.1154962 -2.2459395 -2.3607764 -1.8075186
## theta theta theta theta theta theta theta
## -4.0486716 -2.6845164 -3.0661892 -2.4648441 -3.3381097 -2.7804240 -2.6686100
## ln_kappag
## -5.9132725
# Minimum negative log likelihood
mod$opt$objective
## [1] 2318.892