inla.stack: Data stacking for advanced INLA models

In inbo/INLA: Full Bayesian Analysis of Latent Gaussian Models using Integrated Nested Laplace Approximations

Description

Functions for combining data, effects and observation matrices into inla.stack objects, and extracting information from such objects.

Usage

inla.stack.remove.unused(stack)

inla.stack.compress(stack, remove.unused = TRUE)

inla.stack(..., compress = TRUE, remove.unused = TRUE)

inla.stack.sum(data, A, effects, tag = "", compress = TRUE,
  remove.unused = TRUE)

inla.stack.join(..., compress = TRUE, remove.unused = TRUE)

inla.stack.index(stack, tag)

inla.stack.LHS(stack)

inla.stack.RHS(stack)

inla.stack.data(stack, ...)

inla.stack.A(stack)

Arguments

stack	A inla.data.stack object, created by a call to inla.stack, inla.stack.sum, or inla.stack.join.
remove.unused	If TRUE, compress the model by removing rows of effects corresponding to all-zero columns in the A matrix (and removing those columns).
...	For inla.stack.join, two or more data stacks of class inla.data.stack, created by a call to inla.stack, inla.stack.sum, or inla.stack.join. For inla.stack.data, a list of variables to be joined with the data list.
compress	If TRUE, compress the model by removing duplicated rows of effects, replacing the corresponding A-matrix columns with a single column containing the sum.
data	A list or codedata.frame of named data vectors. Scalars are expanded to match the number of rows in the A matrices, or any non-scalar data vectors. An error is given if the input is inconsistent.
A	A list of observation matrices. Scalars are expanded to diagonal matrices matching the effect vector lengths. An error is given if the input is inconsistent or ambiguous.
effects	A collection of effects/predictors. Each list element corresponds to an observation matrix, and must either be a single vector, a list of vectors, or a data.frame. Single-element effect vectors are expanded to vectors matching the number of columns in the corresponding A matrix. An error is given if the input is inconsistent or ombiguous.
tag	A string specifying a tag for later identification.

Details

For models with a single effects collection, the outer list container for A and effects may be omitted.

Component size definitions:

• [n_l] effect blocks

• [n_k] effects

• [n_i] data values

• [n_jl] effect size for block l

• [n_j] sum_l n_jl total effect size

Input:

• [data] p vectors, each of length n_i

• [A] matrices of size n_i by n_jl

• [effects] collections of effect vectors of length n_jl

predictor(y^1, …, y^p) ~ sum_{l=1}^{n_l} A^l sum_{k=1}^{n_k} g(k, x^{k,l}) = tilde{A} sum_{k=1}^{n_k} g(k, tilde{x}^k)

where

tilde{A} = cbind( A^1, ..., A^{n_l} )

tilde{x}^k = rbind( x^{k,1}, ..., x^{k,n_l} )

and for each block l, any missing x^{k,l} is replaced by an NA vector.

Value

A data stack of class inla.data.stack. Elements:

• data =(y^1, …, y^p, \tilde{x}^1, …, \tilde{x}^{n_k})

• A =\tilde{A}

• data.names List of data names, length p

• effect.names List of effect names, length n_k

• n.data Data length, n_i

• index List indexed by tags, each element indexing into i=1, …, n_i

Functions

• inla.stack.remove.unused: Remove unused entries from an existing stack

• inla.stack.compress: Compress an existing stack by removing duplicates

• inla.stack: Shorthand for inla.stack.join and inla.stack.sum

• inla.stack.sum: Create data stack as a sum of predictors

• inla.stack.join: Join two or more data stacks

• inla.stack.index: Extract tagged indices

• inla.stack.LHS: Extract data associated with the "left hand side" of the model (e.g. the data itself, Ntrials, link, E)

• inla.stack.RHS: Extract data associated with the "right hand side" of the model (all the covariates/predictors)

• inla.stack.data: Extract data for an inla call, and optionally join with other variables

• inla.stack.A: Extract the "A matrix" for control.predictor

See Also

inla.spde.make.A, inla.spde.make.index

Examples

n <- 200
loc <- matrix(runif(n*2), n, 2)
mesh <- inla.mesh.2d(loc.domain = loc,
                     max.edge=c(0.05, 0.2))
proj.obs <- inla.mesh.projector(mesh, loc = loc)
proj.pred <- inla.mesh.projector(mesh, loc = mesh$loc)
spde <- inla.spde2.pcmatern(mesh,
                            prior.range = c(0.01, 0.01),
                            prior.sigma = c(10, 0.01))

covar <- rnorm(n)
field <- inla.qsample(n = 1, Q = inla.spde.precision(spde, theta = log(c(0.5, 1))))[,1]
y <- 2*covar + inla.mesh.project(proj.obs, field)

A.obs <- inla.spde.make.A(mesh, loc = loc)
A.pred = inla.spde.make.A(mesh, loc = proj.pred$loc)
stack.obs <-
    inla.stack(data=list(y=y),
               A=list(A.obs, 1),
               effects=list(c(list(Intercept = 1),
                              inla.spde.make.index("spatial", spde$n.spde)),
                            covar=covar),
               tag="obs")
stack.pred <-
    inla.stack(data=list(y=NA),
               A=list(A.pred),
               effects=list(c(list(Intercept = 1),
                              inla.spde.make.index("spatial", mesh$n))),
               tag="pred")
stack <- inla.stack(stack.obs, stack.pred)

formula <- y ~ -1 + Intercept + covar + f(spatial, model=spde)
result1 <- inla(formula,
                data=inla.stack.data(stack.obs, spde = spde),
                family="gaussian",
                control.predictor = list(A = inla.stack.A(stack.obs),
                                        compute = TRUE))

plot(y, result1$summary.fitted.values[inla.stack.index(stack.obs,"obs")$data, "mean"],
     main = "Observations vs posterior predicted values at the data locations")

result2 <- inla(formula,
                data=inla.stack.data(stack, spde = spde),
                family="gaussian",
                control.predictor = list(A = inla.stack.A(stack),
                                         compute = TRUE))

field.pred <- inla.mesh.project(proj.pred,
  result2$summary.fitted.values[inla.stack.index(stack,"pred")$data, "mean"])
field.pred.sd <- inla.mesh.project(proj.pred,
  result2$summary.fitted.values[inla.stack.index(stack,"pred")$data, "sd"])

plot(field, field.pred, main = "True vs predicted field")
abline(0, 1)
image(inla.mesh.project(mesh,
                        field = field,
                        dims = c(200,200)),
      main = "True field")
image(inla.mesh.project(mesh,
                        field = field.pred,
                        dims = c(200,200)),
      main = "Posterior field mean")
image(inla.mesh.project(mesh,
                        field = field.pred.sd,
                        dims = c(200,200)),
      main = "Prediction standard deviation")
plot(field, (field.pred - field) / 1,
     main = "True field vs standardised prediction residuals")

inbo/INLA documentation built on Dec. 6, 2019, 9:51 a.m.