Beside the basic latent models and likelihoods, INLA provide several ``tools'' to manipulate models and likelihoods. This page describe these.
Very often we encounter the situation where an element of a model is needed more than once for each observation. One example is where
y = a + b*w + ... (*)
for fixed weights 'w' and where (a, b) are bivariate Normals, where (a_i, b_i) is bivariate Normal and all 2-vectors are independent. Using the model
f(idx, model="iid2d", n=2*m, ...)
provide a random vector `v', say, with length 2*m stored as
v = (a_1, a_2, ..., a_m, b_1, b_2, ...., b_m).
To implement the model (*), we simply create an indentical copy of v, v*, where v == v* (nearly). Using the copy-feature, we can spesify the model (*) as
idx = 1:m
idx.copy = m + 1:m
formula = y ~ f(idx, model="iid2d", n=2*m) + f(idx.copy, w, copy="idx") + ....
recalling that the first `m' elements is `a' and the last `m' elements are `b', and where 'w' are the weights.
The second `f()' terms define itself as a copy of `f(idx, ...)', and it inherit some of its features:
n
values
A copied model may also have an unknown scaling (hyperparameter), which is default fixed to be 1. In the following example, we will use this feature to estimate the unknown scaling (in this case, scaling = 2), assuming we know the precision for 'z'.
n=1000
i=1:n
j = i
z = rnorm(n)
w = runif(n)
y = z + 2*z*w + rnorm(n)
formula = y ~ f(i, model="iid",initial=0, fixed=T) +
f(j, w, copy="i", fixed=FALSE)
r = inla(formula, data = data.frame(i,j,w,y))
If the scaling parameter is within given range, then option `range = c(low, high)', can be given. In this case `beta = low + (high-low)*exp(beta.local)/(1+exp(beta.local))', and the prior is defined on beta.local. If low=high or range = NULL, then the identity mapping is used. If high=Inf and low!=Inf, then the mapping `low + exp(beta.local)' is used. The case low=Inf and high!=Inf is not yet implemented.
A model or a copied model can be copied several times. The degree of closeness of 'v' and 'v*' is specified by the argument precision, as the precision of the noise added to 'v' to get 'v*'.
Independent replications of a model with the same hyperparmeters can be defined using the argument 'replicate',
f(idx, model = .., replicate = r)
Here, 'r' is a vector of the same length as 'idx'. In this case, we use a two-dimensional index to index this (sub-)model: (idx, r), so (1,2) identify the first value of the second replication of the model. Number of replications are defined as max(replicate), unless it is defined by the argument 'nrep'.
One example is the model `iid':
i = 1:n
formula = y ~ f(i, model = "iid") + ...
which has an alternative equivalent formulation as `n' replications of an iid-model with length 1
i = rep(1,n)
r = 1:n
formula = y ~ f(i, model="iid", replicate = r) + ...
In the following example, we estimate the parameters in three AR(1) time-series with the same hyperparameters (ie its replicated) but with separate means:
n = 100
y1 = arima.sim(n=n, model=list(ar=c(0.9)))+10
y2 = arima.sim(n=n, model=list(ar=c(0.9)))+20
y3 = arima.sim(n=n, model=list(ar=c(0.9)))+30
formula = y ~ mean -1 + f(i, model="ar1", replicate=r)
y = c(y1,y2,y3)
i = rep(1:n, 3)
r = rep(1:3, each=n)
mean = as.factor(r)
result = inla(formula, family = "gaussian",
data = data.frame(y, i, mean),
control.family = list(initial = 12, fixed=TRUE))
All other arguments is interpreted for the basic model and also replicated. Like argument constr=TRUE, is interpreted as each replication sums to zero, and additional constraints are also replicated.
There is no constraint in INLA that the type of likelihood must be the same for all observations. In fact, every observation could have its own likelihood. Extentions include more than one familily, like the Normal, Poisson, etc, but also having in the model groups of observations with separate hyperparameters within each group, where the family, for example, can be the same.
In the formula
y ~ a + 1
we allow 'y' to be a matrix. In this case each column of `y' define one likelihood where
the family is the same
the hyperparameters are the same.
For each row, only one of the columns could (but don't have to) have an observation (non-NA value), the other colums must have value 'NA'. All other parameters to the likelihood, like
E
Ntrials
offset
scale
are used as appropriate. Example: If row 'i' column 'j' is a Poission observation, then `E[i]' is used as the scaling. Similar with the others. This works as only one column for each row is non-NA.
The argument 'family' is in the case where 'y' is a matrix, a list of families. The argument 'control.family' is then a list of lists; one for each family. OOPS: It is difficult to check that all arguments are internally consistent in this case, therefore strange 'parsing errors' can occur if these arguments are not internally consistent.
To illustrate the ideas, let us do some examples. There are 'real' examples in the case-studies section. (Note that for survival models then the responce 'y' must be a list instead of a matrix, see this example for details.)
The first example, is a simple linear regression, where the first half of the data is observed with unknown precision tau.1 (with a 'default' noninformative prior) and the second half of the data is observed with unknown precision tau.2. Otherwise, the two models have the same form for the linear predictor.
## Simple linear regression with observations with two different
## variances.
n = 100
N = 2*n
y = numeric(N)
x = rnorm(N)
y[1:n] = 1 + x[1:n] + rnorm(n, sd = 1/sqrt(1))
y[1:n + n] = 1 + x[1:n + n] + rnorm(n, sd = 1/sqrt(2))
Y = matrix(NA, N, 2)
Y[1:n, 1] = y[1:n]
Y[1:n + n, 2] = y[1:n + n]
formula = Y ~ x + 1
result = inla(
formula,
data = list(Y=Y, x=x),
family = c("gaussian", "gaussian"),
control.family = list(list(prior = "flat", param = numeric()),
list()))
The second artificial example shows how to use information from two sources to estimate the effect of the covariate 'x'.
## Simple example with two types of likelihoods
n = 10
N = 2*n
## common covariates
x = rnorm(n)
## Poisson, depends on x
E1 = runif(n)
y1 = rpois(n, lambda = E1*exp(x))
## Binomial, depends on x
size = sample(1:10, size=n, replace=TRUE)
prob = exp(x)/(1+exp(x))
y2 = rbinom(n, size= size, prob = prob)
## Join them together
Y = matrix(NA, N, 2)
Y[1:n, 1] = y1
Y[1:n + n, 2] = y2
## The E for the Poisson
E = numeric(N)
E[1:n] = E1
E[1:n + n] = NA
## Ntrials for the Binomial
Ntrials = numeric(N)
Ntrials[1:n] = NA
Ntrials[1:n + n] = size
## Duplicate the covariate which is shared
X = numeric(N)
X[1:n] = x
X[1:n + n] = x
## Formula involving Y as a matrix
formula = Y ~ X - 1
## `family' is now
result = inla(formula,
family = c("poisson", "binomial"),
data = list(Y=Y, X=X),
E = E, Ntrials = Ntrials)
If the covariate 'x' is different for the two families, 'x' and 'xx', say, then we only need to make the following changes
X = numeric(N)
X[1:n] = x
X[1:n + n] = NA
XX = numeric(N)
XX[1:n] = NA
XX[1:n + n] = xx
formula = Y ~ X + XX -1
and add 'XX' into the data.frame. Again this example is artificial as there is no gain of doing a simultanous analysis. However, it shows how the model has to be modified into a 'union' of models with the use of NA's to remove terms.
In the third example, we use also the 'replicate' feature to estimate three replicated AR(1) models with the same hyperparamters, each observed differently.
## An example with three independent AR(1)'s with separate means, but
## with the same hyperparameters. These are observed with three
## different likelihoods.
n = 100
x1 = arima.sim(n=n, model=list(ar=c(0.9))) + 0
x2 = arima.sim(n=n, model=list(ar=c(0.9))) + 1
x3 = arima.sim(n=n, model=list(ar=c(0.9))) + 2
## Binomial observations
Nt = 10 + rpois(n,lambda=1)
y1 = rbinom(n, size=Nt, prob = exp(x1)/(1+exp(x1)))
## Poisson observations
Ep = runif(n, min=1, max=10)
y2 = rpois(n, lambda = Ep*exp(x2))
## Gaussian observations
y3 = rnorm(n, mean=x3, sd=0.1)
## stack these in a 3-column matrix with NA's where not observed
y = matrix(NA, 3*n, 3)
y[1:n, 1] = y1
y[n + 1:n, 2] = y2
y[2*n + 1:n, 3] = y3
## define the model
r = c(rep(1,n), rep(2,n), rep(3,n))
rf = as.factor(r)
i = rep(1:n, 3)
formula = y ~ f(i, model="ar1", replicate=r, constr=TRUE) + rf -1
data = list(y=y, i=i, r=r, rf=rf)
## parameters for the binomial and the poisson
Ntrial = rep(NA, 3*n)
Ntrial[1:n] = Nt
E = rep(NA, 3*n)
E[1:n + n] = Ep
result = inla(formula, family = c("binomial", "poisson", "normal"),
data = data, Ntrial = Ntrial, E = E,
control.family = list(
list(),
list(),
list(initial=0)))
In some cases, the data/response might depend on a linear combination of the 'linear predictor' defined in the formula, like
y ~ 1 + z
then this implies that 'y[1]' depends on 'intercept + beta*z[1]'. Suppose if 'y[1]' depends on '2*intercept + beta*z[1] + beta*z[2]' ? Although it is possible to express this, using the tools we already have, it is more convenient to add another layer into the model. Let 'A' be a 'm x n' matrix, which defines new linear predictors, 'eta~' from 'eta', like
eta~ = A %*% eta
Here, 'eta' is the 'ordinary linear predictor' defined using the formula, and the data depends on the linear predictor 'eta~'. So we might express this as
y ~ 1 + z, with addition matrix A
meaning in short, that
y ~ eta~ (no intercept)
eta~ = A %*% eta
eta = intercept + beta*z
The spesification in inla() is to add the A-matrix using 'control.predictor=list(A=A)'.
The argument 'offset', which might be defined in the formula as 'offset(value)' or as an argument 'inla(..., offset = value)', does have different interpretation in the case where the A-matrix is used. The rule is that 'offset' in the formula, goes into 'eta', whereas an argument 'offset' goes into 'eta~'. If we write out the expressions above adding offsets, 'offset.formula' and 'offset.arg', we get
eta~ = A %*% eta + offset.arg
eta = intercept + beta*z + offset.formula
In the case where there is no A-matrix, then 'offset.total = offset.arg + offset.formula'.
The following example should provide more insight. You may change 'n' and 'm', such that 'm < n', 'm=n' or 'm > n'. Note that since the response has dimension 'm' and the covariates dimension 'n', we need to use 'list(y=y, z=z)' and not a data.frame(). This example also illustrates the use of offset's.
## 'm' is the number of observations of eta*, where eta* = A eta +
## offset.arg, and A is a fixed m x n matrix, and eta has length n. An
## offset in the formula goes into 'eta' whereas an offset in the
## argument of the inla-call, goes into eta*
n = 10
m = 100
offset.formula = 10+ 1:n
offset.arg = 1 + 1:m
## a covariate
z = runif(n)
## the linear predictor eta
eta = 1 + z + offset.formula
## the linear predictor eta* = A eta + offset.arg.
A = matrix(runif(n*m), m, n);
##A = inla.as.sparse(A) ## sparse is ok
## need 'as.vector', as 'Eta' will be a sparseMatrix if 'A' is sparse
## even if ncol(Eta) = 1
Eta = as.vector(A %*% eta) + offset.arg
s = 1e-6
Y = Eta + rnorm(m, sd=s)
## for a check, we can use several offsets. here, m1=-1 and p1=1, so
## they m1+p1 = 0.
r = inla(Y ~ 1+z + offset(offset.formula) + offset(m1) + offset(p1),
## The A-matrix defined here
control.predictor = list(A = A, compute=TRUE, precision = 1e6),
## OOPS: we need to use a list() as the different lengths of Y
## and z
data = list(Y=Y, z=z,
m1 = rep(-1, n),
p1 = rep(1, n),
offset.formula = offset.formula,
offset.arg = offset.arg),
## this is the offset defined in the argument of inla
offset = offset.arg,
##
control.family = list(initial = log(1/s^2), fixed=TRUE))
## this should be a small number
print(max(abs(r$summary.linear.predictor$mean - c(Eta, eta))))
Here is a another example where the informal formula is ``y = intercept + s[j] + 0.5*s[k] + noise''. Instead of using the copy feature, we can implement this model using the A-matrix feature. What we do, is to first define a linear predictor being the intercept and 's', then we use the A-matrix to 'construct the model'. Here is the example:
n = 100L
s = c(-1, 0, 1)
nS = length(s)
j = sample(1L:nS, n, replace=TRUE)
k = j
k[j == 1L] = 2
k[j == 2L] = 3
k[k == 3L] = 1
noise = rnorm(n, sd=0.0001)
y = 1 + s[j] + 0.5*s[k] + noise
## build the formula such that the linear predictor is the intercept
## (index 1) and the 's' term (index 2:(n+1)). then kind of
## 'construct' the model using the A-matrix.
formula = y ~ -1 + intercept + f(idx)
A = matrix(0, n, nS+1L)
for(i in 1L:n) {
A[i, 1L] = 1
A[i, 1L + j[i]] = 1
A[i, 1L + k[i]] = 0.5
}
data = list(intercept = c(1, rep(NA, nS)), idx = c(NA, 1L:nS))
result = inla(formula, data=data, control.predictor=list(A=A))
## should be a straight line
plot(result$summary.random$idx$mean, s)