This document is to record the details of some of the internal specifications of R-INLA.
The purpose of this is largely to document what R-INLA is doing for comparing model implementations in other software.
Make sure to be using a version of R-INLA more recent than INLA_20.06.15 testing version when the internal implementation of the RW1 and RW2 models was updated to be consistent with manual model scaling.
library(INLA)## Loading required package: Matrix## Loading required package: foreach## Loading required package: parallel## Loading required package: sp## This is INLA_22.05.07 built 2022-05-07 09:58:26 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - To enable PARDISO sparse library; see inla.pardiso()Create some very simple space-time data for testing R-INLA.
Data are for four areas in which area 2 is connected to all other areas, areas 1 and 3 are connected, and area 4 is connected only to area 2.
Three time points are simulated.
set.seed(1)
data <- expand.grid(area = 1:4,
time = 1:3)
data$y <- rpois(nrow(data), 2.5)
adj <- rbind(c(0, 1, 1, 0),
c(0, 0, 1, 1),
0,
0)
adj <- adj + t(adj)
rownames(adj) <- colnames(adj) <- letters[1:4]
adj## a b c d
## a 0 1 1 0
## b 1 0 1 1
## c 1 1 0 0
## d 0 1 0 0Structure matrices for ICAR area effect and RW1 time effects.
R_area <- diag(rowSums(adj)) - adj
R_time <- t(diff(diag(3))) %*% diff(diag(3))
R_area## a b c d
## a 2 -1 -1 0
## b -1 3 -1 -1
## c -1 -1 2 0
## d 0 -1 0 1R_time## [,1] [,2] [,3]
## [1,] 1 -1 0
## [2,] -1 2 -1
## [3,] 0 -1 1R-INLA internalsThis section documents details of how certain arguments and options are implemented by R-INLA such as how the small constant is added to the diagonal, model scaling, and Kronecker product for the f(..., group = <>) option.
The strategy for checking is to create a “null” data set with no observations and “fit” the model with fixed values for the hyper parameters to recover the precision matrix constructed by R-INLA.
datanull <- data[data$area == 1, ]
datanull$y <- NA
hyper_area <- list(prec = list(initial = log(1), fixed = TRUE))
hyper_time <- list(prec = list(initial = log(1), fixed = TRUE))fit1 <- inla(y ~ 0 + f(area, model = "besag", graph = adj, hyper = hyper_area),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))INLA:::inla.set.f.default()$diagonal.To see this, in fit1 the fixed value for the precision is 1.0 and in fit2 the value for the precision is 2.0. The added value along the diagonal in both cases is 1e-5.
hyper_area2 <- list(prec = list(initial = log(2), fixed = TRUE))
fit2 <- inla(y ~ 0 + f(area, model = "besag", graph = adj, hyper = hyper_area2),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit1$misc$configs$config[[1]]$Q[-(1:3), -(1:3)]## 4 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 3e+06 -1.0000 -1.0000 .
## [2,] . 3.0001 -1.0000 -1.0000
## [3,] . . 2.0001 .
## [4,] . . . 1.0001fit2$misc$configs$config[[1]]$Q[-(1:3), -(1:3)]## 4 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 3000002 -2.0000 -2.0000 .
## [2,] . 6.0001 -2.0000 -2.0000
## [3,] . . 4.0001 .
## [4,] . . . 2.0001f(..., scale.model = TRUE)When argument scale.model = TRUE, the precision matrix is scaled so that the generalised variance is 1.
fit <- inla(y ~ 0 + f(area, model = "besag", graph = adj, hyper = hyper_area, scale.model = TRUE),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit$misc$configs$constr## $nc
## [1] 1
##
## $A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 0 0 1 1 1 1
##
## $e
## [1] 0fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)]## 4 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 2999999 -0.3565926 -0.3565926 .
## [2,] . 1.0698778 -0.3565926 -0.3565926
## [3,] . . 0.7132852 .
## [4,] . . . 0.3566926inla.scale.model(R_area, constr = list(A = matrix(1, ncol = 4), e = 0))## 4 x 4 sparse Matrix of class "dgTMatrix"
## a b c d
## a 0.7131852 -0.3565926 -0.3565926 .
## b -0.3565926 1.0697778 -0.3565926 -0.3565926
## c -0.3565926 -0.3565926 0.7131852 .
## d . -0.3565926 . 0.3565926Even if constr = FALSE or an alternative constraint is specified in the f() object, the same sum-to-zero constraint is applied to the model scaling.
fit_no_constr <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE, constr = FALSE),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit_no_constr$misc$configs$constr## NULLfit_no_constr$misc$configs$config[[1]]$Q[-(1:3), -(1:3)]## 4 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 2999999 -0.3565926 -0.3565926 .
## [2,] . 1.0697780 -0.3565926 -0.3565926
## [3,] . . 0.7131854 .
## [4,] . . . 0.3565927fit_alt_constr <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE, constr = FALSE,
diagonal = 0,
extraconstr = list(A = matrix(c(1, 1, 0, 0), 1), e = 3)),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit_alt_constr$misc$configs$constr## $nc
## [1] 1
##
## $A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 0 0 1 1 0 0
##
## $e
## [1] 3fit_alt_constr$misc$configs$config[[1]]$Q[-(1:3), -(1:3)]## 4 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 2999999 -0.3565926 -0.3565926 .
## [2,] . 1.0697780 -0.3565926 -0.3565926
## [3,] . . 0.7131854 .
## [4,] . . . 0.3565927Called externally, the alternative constraint does slightly affect model scaling, and so the above confirms that the alternative constraint is not used in the scale.model specification.
inla.scale.model(R_area, constr = list(A = matrix(c(1, 1, 1, 1), ncol = 4), e = 0))## 4 x 4 sparse Matrix of class "dgTMatrix"
## a b c d
## a 0.7131852 -0.3565926 -0.3565926 .
## b -0.3565926 1.0697778 -0.3565926 -0.3565926
## c -0.3565926 -0.3565926 0.7131852 .
## d . -0.3565926 . 0.3565926inla.scale.model(R_area, constr = list(A = matrix(c(1, 1, 0, 0), ncol = 4), e = 3))## 4 x 4 sparse Matrix of class "dgTMatrix"
## a b c d
## a 0.7135650 -0.3567825 -0.3567825 .
## b -0.3567825 1.0703476 -0.3567825 -0.3567825
## c -0.3567825 -0.3567825 0.7135650 .
## d . -0.3567825 . 0.3567825f(..., group = <var>)Next we review the default behaviour of the grouping option to specify product smooths. The below fits a separable space-time model with an ICAR area effect and RW1 time effect, otherwise using defaults.
fit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
group = time, control.group = list(model = "rw1")),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))By default, the a separate sum-to-zero constraint is specified for each level of the group variable.
fit$misc$configs$constr## $nc
## [1] 3
##
## $A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 0 0 0 1 1 1 1 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 1 1 1 1 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 1 1 1
## [,15]
## [1,] 0
## [2,] 0
## [3,] 1
##
## $e
## [1] 0 0 0If extraconstr is specified, the constraint must be the length of the number of levels for the primary variable (e.g. 4 areas in the example below).
This constraint is repeated for each level of the group variable.
area_constr <- list(A = matrix(c(1.5, 0.5, 0, 0), 1), e = 3)
fit_alt_constr <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
group = time, control.group = list(model = "rw1"),
constr = FALSE, extraconstr = area_constr),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit_alt_constr$misc$configs$constr## $nc
## [1] 3
##
## $A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 0 0 0 1.5 0.5 0 0 0.0 0.0 0 0 0.0 0.0 0
## [2,] 0 0 0 0.0 0.0 0 0 1.5 0.5 0 0 0.0 0.0 0
## [3,] 0 0 0 0.0 0.0 0 0 0.0 0.0 0 0 1.5 0.5 0
## [,15]
## [1,] 0
## [2,] 0
## [3,] 0
##
## $e
## [1] 3 3 3As far as I can tell, there is no way to specify (1) a constraint for only some group levels, (2) different constraints for different group levels, or (3) constraints that span multiple group levels. For example, the following constraint with dimension 4 x 3 = 12 might be used to specify an overall sum-to-zero constraint (across all groups) on the space x time latent field.
wishful_constr <- list(A = matrix(1, ncol = 12), e = 0)
wishful_constr## $A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 1 1 1 1 1 1 1 1 1 1 1
##
## $e
## [1] 0But this throws an error because R-INLA is expecting the constraint to have dimension 4 (the number of areas) and will repeat this constraint multiple times.
inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
group = time, control.group = list(model = "rw1"),
constr = FALSE,
extraconstr = wishful_constr),
data = datanull, family = "poisson")## Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
##
## ncol in matrix A(extraconstr) does not correspont to the length of f: 12 4
##
## *** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=1/2
## Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
##
## ncol in matrix A(extraconstr) does not correspont to the length of f: 12 4
##
## *** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=2/2
## Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
##
## ncol in matrix A(extraconstr) does not correspont to the length of f: 12 4## Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts, : *** Fail to get good enough initial values. Maybe it is due to something else.More flexible constraints should be possible by specifying custom model using the "rgeneric" model type and manually specifying the Cmatrix as the Kronecker product.
control.group = list(...) default for scale.model is TRUE, even when the default for scale.model is FALSE.f(..., scale.model = ) default.inla.getOption()$scale.model.default## [1] FALSEModel with defaults for scale model:
fit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
group = time, control.group = list(model = "rw1")),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
R_area_scaled <- inla.scale.model(R_area,
constr = list(A = matrix(1, ncol = 4), e = 0))
R_time_scaled <- inla.scale.model(R_time,
constr = list(A = matrix(1, ncol = 3), e = 0))
round(fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)], 5)## 12 x 12 sparse Matrix of class "dgCMatrix"
##
## [1,] 1e+06 -0.40934 -0.40934 . -0.81867 0.40934 0.40934
## [2,] . 1.22811 -0.40934 -0.40934 0.40934 -1.22801 0.40934
## [3,] . . 0.81877 . 0.40934 0.40934 -0.81867
## [4,] . . . 0.40944 . 0.40934 .
## [5,] . . . . 1000001.07948 -0.81867 -0.81867
## [6,] . . . . . 2.45612 -0.81867
## [7,] . . . . . . 1.63745
## [8,] . . . . . . .
## [9,] . . . . . . .
## [10,] . . . . . . .
## [11,] . . . . . . .
## [12,] . . . . . . .
##
## [1,] . . . . .
## [2,] 0.40934 . . . .
## [3,] . . . . .
## [4,] -0.40934 . . . .
## [5,] . -8.1867e-01 0.40934 0.40934 .
## [6,] -0.81867 4.0934e-01 -1.22801 0.40934 0.40934
## [7,] . 4.0934e-01 0.40934 -0.81867 .
## [8,] 0.81877 . 0.40934 . -0.40934
## [9,] . 1.0000e+06 -0.40934 -0.40934 .
## [10,] . . 1.22811 -0.40934 -0.40934
## [11,] . . . 0.81877 .
## [12,] . . . . 0.40944This matches the Kronecker product of the unscaled area structure matrix and scaled time structure matrix.
kronecker(R_time_scaled, R_area)## 12 x 12 sparse Matrix of class "dgTMatrix"
##
## [1,] 0.8186736 -0.4093368 -0.4093368 . -0.8186736 0.4093368
## [2,] -0.4093368 1.2280105 -0.4093368 -0.4093368 0.4093368 -1.2280105
## [3,] -0.4093368 -0.4093368 0.8186736 . 0.4093368 0.4093368
## [4,] . -0.4093368 . 0.4093368 . 0.4093368
## [5,] -0.8186736 0.4093368 0.4093368 . 1.6373473 -0.8186736
## [6,] 0.4093368 -1.2280105 0.4093368 0.4093368 -0.8186736 2.4560209
## [7,] 0.4093368 0.4093368 -0.8186736 . -0.8186736 -0.8186736
## [8,] . 0.4093368 . -0.4093368 . -0.8186736
## [9,] . . . . -0.8186736 0.4093368
## [10,] . . . . 0.4093368 -1.2280105
## [11,] . . . . 0.4093368 0.4093368
## [12,] . . . . . 0.4093368
##
## [1,] 0.4093368 . . . . .
## [2,] 0.4093368 0.4093368 . . . .
## [3,] -0.8186736 . . . . .
## [4,] . -0.4093368 . . . .
## [5,] -0.8186736 . -0.8186736 0.4093368 0.4093368 .
## [6,] -0.8186736 -0.8186736 0.4093368 -1.2280105 0.4093368 0.4093368
## [7,] 1.6373473 . 0.4093368 0.4093368 -0.8186736 .
## [8,] . 0.8186736 . 0.4093368 . -0.4093368
## [9,] 0.4093368 . 0.8186736 -0.4093368 -0.4093368 .
## [10,] 0.4093368 0.4093368 -0.4093368 1.2280105 -0.4093368 -0.4093368
## [11,] -0.8186736 . -0.4093368 -0.4093368 0.8186736 .
## [12,] . -0.4093368 . -0.4093368 . 0.4093368Fit both scaled f(..., scale.model = TRUE, ..., control.group = list(..., scale.model = TRUE)):
fit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE,
group = time,
control.group = list(model = "rw1", scale.model = TRUE)),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
round(fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)], 5)## 12 x 12 sparse Matrix of class "dgCMatrix"
##
## [1,] 999999.7 -0.14597 -0.14597 . -2.9193e-01 0.14597 0.14597
## [2,] . 0.43800 -0.14597 -0.14597 1.4597e-01 -0.43790 0.14597
## [3,] . . 0.29203 . 1.4597e-01 0.14597 -0.29193
## [4,] . . . 0.14607 . 0.14597 .
## [5,] . . . . 1.0000e+06 -0.29193 -0.29193
## [6,] . . . . . 0.87590 -0.29193
## [7,] . . . . . . 0.58397
## [8,] . . . . . . .
## [9,] . . . . . . .
## [10,] . . . . . . .
## [11,] . . . . . . .
## [12,] . . . . . . .
##
## [1,] . . . . .
## [2,] 0.14597 . . . .
## [3,] . . . . .
## [4,] -0.14597 . . . .
## [5,] . -0.29193 0.14597 0.14597 .
## [6,] -0.29193 0.14597 -0.43790 0.14597 0.14597
## [7,] . 0.14597 0.14597 -0.29193 .
## [8,] 0.29203 . 0.14597 . -0.14597
## [9,] . 999999.73407 -0.14597 -0.14597 .
## [10,] . . 0.43800 -0.14597 -0.14597
## [11,] . . . 0.29203 .
## [12,] . . . . 0.14607round(kronecker(R_time_scaled, R_area_scaled), 5)## 12 x 12 sparse Matrix of class "dgTMatrix"
##
## [1,] 0.29193 -0.14597 -0.14597 . -0.29193 0.14597 0.14597 .
## [2,] -0.14597 0.43790 -0.14597 -0.14597 0.14597 -0.43790 0.14597 0.14597
## [3,] -0.14597 -0.14597 0.29193 . 0.14597 0.14597 -0.29193 .
## [4,] . -0.14597 . 0.14597 . 0.14597 . -0.14597
## [5,] -0.29193 0.14597 0.14597 . 0.58387 -0.29193 -0.29193 .
## [6,] 0.14597 -0.43790 0.14597 0.14597 -0.29193 0.87580 -0.29193 -0.29193
## [7,] 0.14597 0.14597 -0.29193 . -0.29193 -0.29193 0.58387 .
## [8,] . 0.14597 . -0.14597 . -0.29193 . 0.29193
## [9,] . . . . -0.29193 0.14597 0.14597 .
## [10,] . . . . 0.14597 -0.43790 0.14597 0.14597
## [11,] . . . . 0.14597 0.14597 -0.29193 .
## [12,] . . . . . 0.14597 . -0.14597
##
## [1,] . . . .
## [2,] . . . .
## [3,] . . . .
## [4,] . . . .
## [5,] -0.29193 0.14597 0.14597 .
## [6,] 0.14597 -0.43790 0.14597 0.14597
## [7,] 0.14597 0.14597 -0.29193 .
## [8,] . 0.14597 . -0.14597
## [9,] 0.29193 -0.14597 -0.14597 .
## [10,] -0.14597 0.43790 -0.14597 -0.14597
## [11,] -0.14597 -0.14597 0.29193 .
## [12,] . -0.14597 . 0.14597Fit main effect scaled and group unscaled f(..., scale.model = TRUE, ..., control.group = list(..., scale.model = FALSE)).
fit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE,
group = time,
control.group = list(model = "rw1", scale.model = FALSE)),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
round(fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)], 5)## 12 x 12 sparse Matrix of class "dgCMatrix"
##
## [1,] 1e+06 -0.35659 -0.35659 . -0.71319 0.35659 0.35659
## [2,] . 1.06988 -0.35659 -0.35659 0.35659 -1.06978 0.35659
## [3,] . . 0.71329 . 0.35659 0.35659 -0.71319
## [4,] . . . 0.35669 . 0.35659 .
## [5,] . . . . 1000000.86851 -0.71319 -0.71319
## [6,] . . . . . 2.13966 -0.71319
## [7,] . . . . . . 1.42647
## [8,] . . . . . . .
## [9,] . . . . . . .
## [10,] . . . . . . .
## [11,] . . . . . . .
## [12,] . . . . . . .
##
## [1,] . . . . .
## [2,] 0.35659 . . . .
## [3,] . . . . .
## [4,] -0.35659 . . . .
## [5,] . -7.1319e-01 0.35659 0.35659 .
## [6,] -0.71319 3.5659e-01 -1.06978 0.35659 0.35659
## [7,] . 3.5659e-01 0.35659 -0.71319 .
## [8,] 0.71329 . 0.35659 . -0.35659
## [9,] . 1.0000e+06 -0.35659 -0.35659 .
## [10,] . . 1.06988 -0.35659 -0.35659
## [11,] . . . 0.71329 .
## [12,] . . . . 0.35669round(kronecker(R_time, R_area_scaled), 5)## 12 x 12 sparse Matrix of class "dgTMatrix"
##
## [1,] 0.71319 -0.35659 -0.35659 . -0.71319 0.35659 0.35659 .
## [2,] -0.35659 1.06978 -0.35659 -0.35659 0.35659 -1.06978 0.35659 0.35659
## [3,] -0.35659 -0.35659 0.71319 . 0.35659 0.35659 -0.71319 .
## [4,] . -0.35659 . 0.35659 . 0.35659 . -0.35659
## [5,] -0.71319 0.35659 0.35659 . 1.42637 -0.71319 -0.71319 .
## [6,] 0.35659 -1.06978 0.35659 0.35659 -0.71319 2.13956 -0.71319 -0.71319
## [7,] 0.35659 0.35659 -0.71319 . -0.71319 -0.71319 1.42637 .
## [8,] . 0.35659 . -0.35659 . -0.71319 . 0.71319
## [9,] . . . . -0.71319 0.35659 0.35659 .
## [10,] . . . . 0.35659 -1.06978 0.35659 0.35659
## [11,] . . . . 0.35659 0.35659 -0.71319 .
## [12,] . . . . . 0.35659 . -0.35659
##
## [1,] . . . .
## [2,] . . . .
## [3,] . . . .
## [4,] . . . .
## [5,] -0.71319 0.35659 0.35659 .
## [6,] 0.35659 -1.06978 0.35659 0.35659
## [7,] 0.35659 0.35659 -0.71319 .
## [8,] . 0.35659 . -0.35659
## [9,] 0.71319 -0.35659 -0.35659 .
## [10,] -0.35659 1.06978 -0.35659 -0.35659
## [11,] -0.35659 -0.35659 0.71319 .
## [12,] . -0.35659 . 0.35659Cmatrix.fit <- inla(y ~ 1 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE,
group = time,
control.group = list(model = "rw1", scale.model = FALSE)),
data = data, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.compute = list(config = TRUE))
Q <- kronecker(R_time, R_area_scaled)
data$area.time <- 1:12
A <- rbind(rep(c(1, 0), c(4, 8)),
rep(c(0, 1, 0), c(4, 4, 4)),
rep(c(0, 1), c(8, 4)))
e <- c(0, 0, 0)
diagval <- INLA:::inla.set.f.default()$diagonal
fitQ <- inla(y ~ 1 +
f(area.time, model = "generic0", Cmatrix = Q, hyper = hyper_area,
diagonal = diagval, extraconstr = list(A = A, e = e)),
data = data, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.compute = list(config = TRUE))The effective parameters are equal:
fit$neffp## NULLfitQ$neffp## NULLFixed effect and random effects are equal:
fit$summary.fixed[ , 1:2]fitQ$summary.fixed[ , 1:2]fit$summary.random[[1]][ , 1:3]fitQ$summary.random[[1]][ , 1:3]Marginal likelihood slightly different. Suspect difference in scaling constants?
fit$mlik## [,1]
## log marginal-likelihood (integration) -23.1046
## log marginal-likelihood (Gaussian) -23.1046fitQ$mlik## [,1]
## log marginal-likelihood (integration) -25.86141
## log marginal-likelihood (Gaussian) -25.86141The raw AR1 structure matrix is scaled by \(1/(1-\rho^2)\) such that the marginal variance of the precision matrix is the specified precision.
rho <- 0.8The AR1 precision matrix is given by:
R_ar1 <- sparseMatrix(1:5, 1:5, x = 1)
diag(R_ar1)[2:4] <- 1 + rho^2
R_ar1[cbind(1:4, 2:5)] <- -rho
R_ar1[cbind(2:5, 1:4)] <- -rhoFor a marginal precision of 1.0, scale the raw structure matrix:
Q_ar1 <- R_ar1 * 1 / (1 - rho^2)Show this matches the Q matrix constructed by R-INLA:
hyper_ar1 <- list(prec = list(initial = log(1), fixed = TRUE),
rho = list(initial = log((1 + rho) / (1 - rho)), fixed = TRUE))
fit <- inla(y ~ 0 +
f(time, model = "ar1", values = 1:5, hyper = hyper_ar1),
data = datanull[1, ], family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
fit$misc$configs$config[[1]]$Q[-1, -1]## 5 x 5 sparse Matrix of class "dgCMatrix"
##
## [1,] 1000002 -2.222223 . . .
## [2,] . 4.555557 -2.222223 . .
## [3,] . . 4.555557 -2.222223 .
## [4,] . . . 4.555557 -2.222223
## [5,] . . . . 2.777779Q_ar1## 5 x 5 sparse Matrix of class "dgCMatrix"
##
## [1,] 2.777778 -2.222222 . . .
## [2,] -2.222222 4.555556 -2.222222 . .
## [3,] . -2.222222 4.555556 -2.222222 .
## [4,] . . -2.222222 4.555556 -2.222222
## [5,] . . . -2.222222 2.777778Three time points:
R_ar1 <- sparseMatrix(1:3, 1:3, x = 1)
diag(R_ar1)[2] <- 1 + rho^2
R_ar1[cbind(1:2, 2:3)] <- -rho
R_ar1[cbind(2:3, 1:2)] <- -rho
Q_ar1 <- R_ar1 * 1 / (1 - rho^2)
fit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
group = time,
control.group = list(model = "ar1", hyper = hyper_ar1["rho"])),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
round(fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)], 5)## 12 x 12 sparse Matrix of class "dgCMatrix"
##
## [1,] 1000005 -2.77778 -2.77778 . -4.44445 2.22222 2.22222
## [2,] . 8.33344 -2.77778 -2.77778 2.22222 -6.66667 2.22222
## [3,] . . 5.55566 . 2.22222 2.22222 -4.44445
## [4,] . . . 2.77788 . 2.22222 .
## [5,] . . . . 1000008.55325 -4.55556 -4.55556
## [6,] . . . . . 13.66677 -4.55556
## [7,] . . . . . . 9.11121
## [8,] . . . . . . .
## [9,] . . . . . . .
## [10,] . . . . . . .
## [11,] . . . . . . .
## [12,] . . . . . . .
##
## [1,] . . . . .
## [2,] 2.22222 . . . .
## [3,] . . . . .
## [4,] -2.22222 . . . .
## [5,] . -4.44445 2.22222 2.22222 .
## [6,] -4.55556 2.22222 -6.66667 2.22222 2.22222
## [7,] . 2.22222 2.22222 -4.44445 .
## [8,] 4.55566 . 2.22222 . -2.22222
## [9,] . 1000004.99769 -2.77778 -2.77778 .
## [10,] . . 8.33344 -2.77778 -2.77778
## [11,] . . . 5.55566 .
## [12,] . . . . 2.77788kronecker(Q_ar1, R_area)## 12 x 12 sparse Matrix of class "dgTMatrix"
##
## [1,] 5.555556 -2.777778 -2.777778 . -4.444444 2.222222 2.222222
## [2,] -2.777778 8.333333 -2.777778 -2.777778 2.222222 -6.666667 2.222222
## [3,] -2.777778 -2.777778 5.555556 . 2.222222 2.222222 -4.444444
## [4,] . -2.777778 . 2.777778 . 2.222222 .
## [5,] -4.444444 2.222222 2.222222 . 9.111111 -4.555556 -4.555556
## [6,] 2.222222 -6.666667 2.222222 2.222222 -4.555556 13.666667 -4.555556
## [7,] 2.222222 2.222222 -4.444444 . -4.555556 -4.555556 9.111111
## [8,] . 2.222222 . -2.222222 . -4.555556 .
## [9,] . . . . -4.444444 2.222222 2.222222
## [10,] . . . . 2.222222 -6.666667 2.222222
## [11,] . . . . 2.222222 2.222222 -4.444444
## [12,] . . . . . 2.222222 .
##
## [1,] . . . . .
## [2,] 2.222222 . . . .
## [3,] . . . . .
## [4,] -2.222222 . . . .
## [5,] . -4.444444 2.222222 2.222222 .
## [6,] -4.555556 2.222222 -6.666667 2.222222 2.222222
## [7,] . 2.222222 2.222222 -4.444444 .
## [8,] 4.555556 . 2.222222 . -2.222222
## [9,] . 5.555556 -2.777778 -2.777778 .
## [10,] 2.222222 -2.777778 8.333333 -2.777778 -2.777778
## [11,] . -2.777778 -2.777778 5.555556 .
## [12,] -2.222222 . -2.777778 . 2.777778scale.model = TRUEfit <- inla(y ~ 0 +
f(area, model = "besag", graph = adj, hyper = hyper_area,
scale.model = TRUE,
group = time,
control.group = list(model = "ar1", hyper = hyper_ar1["rho"])),
data = datanull, family = "poisson",
control.inla = list(strategy = "gaussian", int.strategy = "eb"),
control.fixed = list(mean.intercept = 0, prec.intercept = 1),
control.compute = list(config = TRUE))
round(fit$misc$configs$config[[1]]$Q[-(1:3), -(1:3)], 5)## 12 x 12 sparse Matrix of class "dgCMatrix"
##
## [1,] 1000001 -0.99054 -0.99054 . -1.58486 0.79243 0.79243
## [2,] . 2.97171 -0.99054 -0.99054 0.79243 -2.37729 0.79243
## [3,] . . 1.98117 . 0.79243 0.79243 -1.58486
## [4,] . . . 0.99064 . 0.79243 .
## [5,] . . . . 1000002.69109 -1.62448 -1.62448
## [6,] . . . . . 4.87353 -1.62448
## [7,] . . . . . . 3.24906
## [8,] . . . . . . .
## [9,] . . . . . . .
## [10,] . . . . . . .
## [11,] . . . . . . .
## [12,] . . . . . . .
##
## [1,] . . . . .
## [2,] 0.79243 . . . .
## [3,] . . . . .
## [4,] -0.79243 . . . .
## [5,] . -1.58486 0.79243 0.79243 .
## [6,] -1.62448 0.79243 -2.37729 0.79243 0.79243
## [7,] . 0.79243 0.79243 -1.58486 .
## [8,] 1.62458 . 0.79243 . -0.79243
## [9,] . 1000001.42321 -0.99054 -0.99054 .
## [10,] . . 2.97171 -0.99054 -0.99054
## [11,] . . . 1.98117 .
## [12,] . . . . 0.99064round(kronecker(Q_ar1, R_area_scaled), 5)## 12 x 12 sparse Matrix of class "dgTMatrix"
##
## [1,] 1.98107 -0.99053 -0.99053 . -1.58486 0.79243 0.79243 .
## [2,] -0.99053 2.97160 -0.99053 -0.99053 0.79243 -2.37728 0.79243 0.79243
## [3,] -0.99053 -0.99053 1.98107 . 0.79243 0.79243 -1.58486 .
## [4,] . -0.99053 . 0.99053 . 0.79243 . -0.79243
## [5,] -1.58486 0.79243 0.79243 . 3.24895 -1.62448 -1.62448 .
## [6,] 0.79243 -2.37728 0.79243 0.79243 -1.62448 4.87343 -1.62448 -1.62448
## [7,] 0.79243 0.79243 -1.58486 . -1.62448 -1.62448 3.24895 .
## [8,] . 0.79243 . -0.79243 . -1.62448 . 1.62448
## [9,] . . . . -1.58486 0.79243 0.79243 .
## [10,] . . . . 0.79243 -2.37728 0.79243 0.79243
## [11,] . . . . 0.79243 0.79243 -1.58486 .
## [12,] . . . . . 0.79243 . -0.79243
##
## [1,] . . . .
## [2,] . . . .
## [3,] . . . .
## [4,] . . . .
## [5,] -1.58486 0.79243 0.79243 .
## [6,] 0.79243 -2.37728 0.79243 0.79243
## [7,] 0.79243 0.79243 -1.58486 .
## [8,] . 0.79243 . -0.79243
## [9,] 1.98107 -0.99053 -0.99053 .
## [10,] -0.99053 2.97160 -0.99053 -0.99053
## [11,] -0.99053 -0.99053 1.98107 .
## [12,] . -0.99053 . 0.99053sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] INLA_22.05.07 sp_1.4-7 foreach_1.5.2 Matrix_1.4-1
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.13 knitr_1.39 magrittr_2.0.3 splines_4.2.0
## [5] lattice_0.20-45 R6_2.5.1 rlang_1.0.2 fastmap_1.1.0
## [9] stringr_1.4.0 tools_4.2.0 grid_4.2.0 xfun_0.31
## [13] cli_3.3.0 jquerylib_0.1.4 MatrixModels_0.5-0 htmltools_0.5.2
## [17] iterators_1.0.14 yaml_2.3.5 digest_0.6.29 bookdown_0.26
## [21] vctrs_0.4.1 sass_0.4.1 codetools_0.2-18 evaluate_0.15
## [25] rmarkdown_2.14 stringi_1.7.6 compiler_4.2.0 bslib_0.3.1
## [29] jsonlite_1.8.0