inla.stack: Data stacking for advanced INLA models

Description Usage Arguments Details Value See Also Examples

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
## 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

Input:

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:

See Also

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

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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.