inla.spde2.pcmatern {INLA} | R Documentation

## Matern SPDE model object with PC prior for INLA

### Description

Create an `inla.spde2` model object for a Matern model, using a PC prior for the parameters.

### Usage

    inla.spde2.pcmatern(mesh,
                        alpha = 2,
                        param = NULL,
                        constr = FALSE,
                        extraconstr.int = NULL,
                        extraconstr = NULL,
                        fractional.method = c("parsimonious", "null"),
                        n.iid.group = 1,
                        prior.range = NULL,
                        prior.sigma = NULL)

### Arguments

`mesh` | The mesh to build the model on, as an `inla.mesh` or `inla.mesh.1d` object.
`alpha` | Fractional operator order, 0 < alpha <= 2 supported, for ν=α-d/2>0.
`param` | Further model parameters. Not currently used.
`constr` | If `TRUE`, apply an integrate-to-zero constraint. Default `FALSE`.
`extraconstr.int` | Field integral constraints.
`extraconstr` | Direct linear combination constraints on the basis weights.
`fractional.method` | Specifies the approximation method to use for fractional (non-integer) `alpha` values. `'parsimonious'` gives an overall approximate minimal covariance error, `'null'` uses approximates low-order properties.
`n.iid.group` | If greater than 1, build an explicitly iid replicated model, to support constraints applied to the combined replicates, for example in a time-replicated spatial model. Constraints can either be specified for a single mesh, in which case it's applied to the average of the replicates (`ncol(A)` should be `mesh$n` for 2D meshes, `mesh$m` for 1D), or as general constraints on the collection of replicates (`ncol(A)` should be `mesh$n * n.iid.group` for 2D meshes, `mesh$m * n.iid.group` for 1D).
`prior.range` | A length 2 vector, with `(range0,Prange)` specifying that P(ρ < ρ_0)=p_ρ, where ρ is the spatial range of the random field. If `Prange` is `NA`, then `range0` is used as a fixed range value.
`prior.sigma` | A length 2 vector, with `(sigma0,Psigma)` specifying that P(σ > σ_0)=p_σ, where σ is the marginal standard deviation of the field. If `Psigma` is `NA`, then `sigma0` is used as a fixed range value.

### Details

This method constructs a Matern SPDE model, with spatial range ρ and standard deviation parameter σ. In the parameterisation

(kappa^2-Delta)^(alpha/2) (tau x(u)) = W(u)

the spatial scale parameter κ=√{8ν}/ρ, where ν=α-d/2, and τ is proportional to 1/σ.

Stationary models are supported for 0 < alpha <= 2, with spectral approximation methods used for non-integer α, with approximation method determined by `fractional.method`.

Integration and other general linear constraints are supported via the `constr`, `extraconstr.int`, and `extraconstr` parameters, which also interact with `n.iid.group`.

The joint PC prior density for the spatial range, ρ, and the marginal standard deviation, σ, and is

p(rho, sigma) = (d R)/2 rho^(-1-d/2) exp(-R rho^(-d/2)) S exp(-S sigma)

where R and S are hyperparameters that must be determined by the analyst. The practical approach for this in INLA is to require the user to indirectly specify these hyperparameters through

P(ρ < ρ_0) = p_ρ

and

P(σ > σ_0) = p_σ

where the user specifies the lower tail quantile and probability for the range (ρ_0 and p_ρ) and the upper tail quantile and probability for the standard deviation (σ_0 and α_σ).

This allows the user to control the priors of the parameters by supplying knowledge of the scale of the problem. What is a reasonable upper magnitude for the spatial effect and what is a reasonable lower scale at which the spatial effect can operate? The shape of the prior was derived through a construction that shrinks the spatial effect towards a base model of no spatial effect in the sense of distance measured by Kullback-Leibler divergence.

The prior is constructed in two steps, under the idea that having a spatial field is an extension of not having a spatial field. First, a spatially constant random effect (ρ = ∞) with finite variance is more complex than not having a random effect (σ = 0). Second, a spatial field with spatial variation (ρ < ∞) is more complex than the random effect with no spatial variation. Each of these extensions are shrunk towards the simpler model and, as a result, we shrink the spatial field towards the base model of no spatial variation and zero variance (ρ = ∞ and σ = 0).

The details behind the construction of the prior is presented in Fuglstad, et al. (2016) and is based on the PC prior framework (Simpson, et al., 2015).

### Value

An `inla.spde2` object.

### Author(s)

Finn Lindgren finn.lindgren@gmail.com

### References

Fuglstad, G.-A., Simpson, D., Lindgren, F., and Rue, H. (2016) Constructing Priors that Penalize the Complexity of Gaussian Random Fields. arXiv:1503.00256

Simpson, D., Rue, H., Martins, T., Riebler, A., and Sørbye, S. (2015) Penalising model component complexity: A principled, practical approach to constructing priors. arXiv:1403.4630

### See Also

`inla.mesh.2d`, `inla.mesh.create`, `inla.mesh.1d`, `inla.mesh.basis`, `inla.spde2.matern`, `inla.spde2.generic`

### Examples

    ## Spatial interpolation
      n = 100
      field.fcn = function(loc) (10*cos(2*pi*2*(loc[,1]+loc[,2])))
      loc = matrix(runif(n*2),n,2)
      ## One field, 2 observations per location
      idx.y = rep(1:n,2)
      y = field.fcn(loc[idx.y,]) + rnorm(length(idx.y))

      mesh = inla.mesh.2d(loc, max.edge=0.05, cutoff=0.01)
      spde = inla.spde2.pcmatern(mesh,
               prior.range=c(0.01,0.1), prior.sigma=c(100,0.1))
      data = list(y=y, field=mesh$idx$loc[idx.y])
      formula = y ~ -1 + f(field, model=spde)
      result = inla(formula, data=data, family="normal")

      ## Plot the mesh structure:
      plot(mesh)

      if (require(rgl)) {
        col.pal = colorRampPalette(c("blue","cyan","green","yellow","red"))
        ## Plot the posterior mean:
        plot(mesh, rgl=TRUE,
             result$summary.random$field[,"mean"],
             color.palette = col.pal)
        ## Plot residual field:
        plot(mesh, rgl=TRUE,
             result$summary.random$field[,"mean"]-field.fcn(mesh$loc),
             color.palette = col.pal)
      }

      result.field = inla.spde.result(result, "field", spde)
      par(mfrow=c(2,1))
      plot(result.field$marginals.range.nominal[[1]],
           type="l", main="Posterior density for range")
      plot(inla.tmarginal(sqrt, result.field$marginals.variance.nominal[[1]]),
           type="l", main="Posterior density for std.dev.")
      par(mfrow=c(1,1))

    ## Spatial model
      set.seed(1234234)

      ## Generate spatial locations
      nObs = 200
      loc = matrix(runif(nObs*2), nrow = nObs, ncol = 2)

      ## Generate observation of spatial field
      nu = 1.0
      rhoT = 0.2
      kappaT = sqrt(8*nu)/rhoT
      sigT = 1.0
      Sig = sigT^2*inla.matern.cov(nu = nu,
                                   kappa = kappaT,
                                   x = as.matrix(dist(loc)),
                                   d = 2,
                                   corr = TRUE)
      L = t(chol(Sig))
      u = L %*% rnorm(nObs)

      ## Construct observation with nugget
      sigN = 0.1
      y = u + sigN*rnorm(nObs)

      ## Create the mesh and spde object
      mesh = inla.mesh.2d(loc,
                          max.edge = 0.05,
                          cutoff = 0.01)
      spde = inla.spde2.pcmatern(mesh,
                                 prior.range = c(0.01, 0.05),
                                 prior.sigma = c(10, 0.05))

      ## Create projection matrix for observations
      A = inla.spde.make.A(mesh = mesh,
                           loc = loc)

      ## Run model without any covariates
      idx = 1:spde$n.spde
      res = inla(y ~ f(idx, model = spde) - 1,
                 data = list(y = y, idx = idx, spde = spde),
                 control.predictor = list(A = A))

    ## Re-run model with fixed range
      spde.fixed = inla.spde2.pcmatern(mesh,
                                       prior.range = c(0.2, NA),
                                       prior.sigma = c(10, 0.05))

      res.fixed = inla(y ~ f(idx, model = spde) - 1,
                       data = list(y = y, idx = idx, spde = spde.fixed),
                       control.predictor = list(A = A))

[Package INLA version 18.07.12 Index]