inla.stack: Data stacking for advanced INLA models

In andrewzm/INLA: Functions which allow to perform full Bayesian analysis of latent Gaussian models using Integrated Nested Laplace Approximaxion

Description

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

Usage

## Create data stack as a sum of predictors:
inla.stack.sum(data, A, effects, tag="", compress=TRUE, remove.unused=TRUE)

## Join two or more data stacks:
inla.stack.join(..., compress=TRUE, remove.unused=TRUE)

## Shorthand for inla.stack.join and inla.stack.sum:
inla.stack(..., compress=TRUE, remove.unused=TRUE)

## Compress an existing stack:
inla.stack.compress(stack, remove.unused=TRUE)

## Remove unused entries from an existing stack:
inla.stack.remove.unused(stack)

## Extract tagged indices:
inla.stack.index(stack, tag)

## Extract data for inla call, and optionally join with other variables:
inla.stack.data(stack, ...)

## Extract "A matrix" for control.predictor:
inla.stack.A(stack)

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

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

Arguments

data	A list of data vectors.
A	A list of observation matrices.
effects	A collection of effects/predictors. Each list element corresponds to an observation matrix, and must either be a single vector or a list of vectors.
tag	A string specifying a tag for later identification.
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.
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).
stack	An inla.data.stack object, created by a call to inla.stack, inla.stack.sum, or inla.stack.join.
...	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.

Details

• [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

See Also

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

Examples

n = 200
loc = matrix(runif(n*2),n,2)
mesh = inla.mesh.create.helper(points.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.matern(mesh,
                         B.tau=cbind(log(1), 1, 0),
                         B.kappa=matrix(c(log(sqrt(8)/0.2), 0, 1), 1, 3))

covar = rnorm(n)
field = inla.qsample(n=1, Q=inla.spde2.precision(spde, theta=c(0,0)))[,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=rep(1, mesh$n)),
                              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=rep(1, mesh$n)),
                              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),
               verbose=TRUE)

plot(y,result1$summary.fitted.values[inla.stack.index(stack.obs,"obs")$data,"mean"])

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

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

dev.new()
plot(field, field.pred)
dev.new()
image(inla.mesh.project(mesh,
    field=field,
    dims=c(200,200)))
dev.new()
image(inla.mesh.project(mesh,
    field=field.pred,
    dims=c(200,200)))

andrewzm/INLA documentation built on May 10, 2019, 11:12 a.m.