Description Usage Arguments Details Value Prior distributions of the hyperparameters References Examples
Multivariate generalization of the intrinsic conditional autorregresive model. The matrix which models the variability between diseases is a symmetric matrix with the inverse of the marginal precisions on the diagonal elements and the correlation parameters divided by the square root of the precisions on the off-diagonal elements.
1 2 3 | inla.rgeneric.IMCAR.model(cmd, theta)
inla.IMCAR.model(...)
|
... |
Arguments to be passed to 'inla.rgeneric.define'. |
cmd |
Arguments used by latent effects defined using the 'rgeneric' latent effect. |
theta |
Vector of hyperparameters. |
This function is used to define a latent effect that is a multivariate spatial effect with a intrinsic conditional autorregresive distribution and a symmetric matrix in order to model the whitin-disease and the between-diseases variability, respectively. Due to this effect is a multivariate spatial latent effect this function requires the following arguments when defining the latent effect:
W Adjacency SPARSE matrix for spatial effect in the basic binary code.
k Number of diseases of the multivariate study.
This model is defined using the 'f()' function and an index in order to identify the spatial areas. See the example.
This is used internally by the 'INLA::inla()'.
The hyperparamenters of this lattent effect are the marginal precisions of each disease which are equal to the number of diseases and the correlation parameters for the whole pair of diseases.
Palmí-Perales F, Gómez-Rubio V, Martinez-Beneito MA (2021). “Bayesian Multivariate Spatial Models for Lattice Data with INLA.” _Journal of Statistical Software_, *98*(2), 1-29. doi: 10.18637/jss.v098.i02 (URL: https://doi.org/10.18637/jss.v098.i02).
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 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 | if (require("INLA", quietly = TRUE)) {
require(spdep)
require(spData)
require(rgdal)
#Load SIDS data
nc.sids <- readOGR(system.file("shapes/sids.shp", package="spData")[1])
proj4string(nc.sids) <- CRS("+proj=longlat +ellps=clrk66")
#Compute adjacency matrix, as nb object 'adj' and sparse matrix 'W'
adj <- poly2nb(nc.sids)
W <- as(nb2mat(adj, style = "B"), "Matrix")
#Compute expected cases
r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
nc.sids$EXP74 <- r74 * nc.sids$BIR74
nc.sids$SMR74 <- nc.sids$SID74 / nc.sids$EXP74
nc.sids$NWPROP74 <- nc.sids$NWBIR74 / nc.sids$BIR74
r79 <- sum(nc.sids$SID79) / sum(nc.sids$BIR79)
nc.sids$EXP79 <- r79 * nc.sids$BIR79
nc.sids$SMR79 <- nc.sids$SID79 / nc.sids$EXP79
nc.sids$NWPROP79 <- nc.sids$NWBIR79 / nc.sids$BIR79
# Define sum-to-zero constraints
A <- kronecker(Diagonal(2, 1), Matrix(1, ncol = nrow(W), nrow = 1))
e = rep(0, 2)
# Data (replicated to assess scalability)
#Real data
n.rep <- 1
d <- list(OBS = c(nc.sids$SID74, nc.sids$SID79),
NWPROP = c(nc.sids$NWPROP74, nc.sids$NWPROP79),
EXP = c(nc.sids$EXP74, nc.sids$EXP79))
d <- lapply(d, function(X) { rep(X, n.rep)})
d$idx <- 1:length(d$OBS)
# Model parameters
k <- 2 * n.rep #Number of diseases
#Define model IMCAR
model <- inla.rgeneric.define(inla.rgeneric.IMCAR.model, debug = FALSE,
k = k, W = W)
#Fit model
r <- inla(OBS ~ 1 + f(idx, model = model,
extraconstr = list(A = as.matrix(A), e = e)), # + NWPROP,
data = d, E = EXP, family = "poisson",
control.compute = list(config = TRUE),
control.predictor = list(compute = TRUE))
summary(r)
# Transformed parameters
r.hyperpar <- inla.MCAR.transform(r, k = 2, model = "IMCAR")
r.hyperpar$summary.hyperpar
#Get fitted data, i.e., relative risk
nc.sids$FITTED74 <- r$summary.fitted.values[1:100, "mean"]
nc.sids$FITTED79 <- r$summary.fitted.values[100 + 1:100, "mean"]
#Display fitted relative risks
dev.new()
spplot(nc.sids, c("SMR74", "FITTED74", "SMR79", "FITTED79"))
#Show marginals of tau_1, tau_2, rho
marg.tau1 <- inla.tmarginal(
function(x) exp(x),
r$marginals.hyperpar[[1]])
marg.tau2 <- inla.tmarginal(
function(x) exp(x),
r$marginals.hyperpar[[2]])
marg.rho <- inla.tmarginal(
function(x) (2*exp(x))/(1 + exp(x)) - 1,
r$marginals.hyperpar[[3]])
dev.new()
oldpar <- par(mfrow = c(2, 2))
plot(marg.tau1, main = "tau1", type = "l")
plot(marg.tau2, main = "tau2", type = "l")
plot(marg.rho, main = "rho", type = "l")
par(oldpar)
## Running UNIVARIATE MODEL
#Real data
n.rep <- 1
d <- list(OBS = nc.sids$SID74,
NWPROP = nc.sids$NWPROP74,
EXP = nc.sids$EXP74)
d <- lapply(d, function(X) { rep(X, n.rep)})
d$idx <- 1:length(d$OBS)
#Fit model
r.uni <- inla(OBS ~ 1 + f(idx, model = "besag", graph = W), # + NWPROP,
data = d, E = EXP, family = "poisson",
control.predictor = list(compute = TRUE))
summary(r.uni)
nc.sids$FITTED74.uni <- r.uni$summary.fitted.values[ , "mean"]
#Display univariate VS multivariate fitted relative risks.
dev.new()
spplot(nc.sids, c("SMR74", "FITTED74", "FITTED74.uni"))
spplot(nc.sids, c("FITTED74", "FITTED74.uni"),
main=list(label="Relative risk estimation",cex=2))
dev.new()
plot(nc.sids$FITTED74.uni, nc.sids$FITTED74,
main="Relative Risk estimations", xlab="Univariate RR estimations",
ylab="Multivariate RR estimations", xlim=c(0.5, 2.5), ylim=c(0.5, 2.5))
abline(h=0, col="grey")
abline(v=0, col="grey")
abline(a=0, b=1, col="red")
#Plot posterior mean of the spatial effects univ VS multi
nc.sids$m.uni <- r.uni$summary.random$idx[, "mean"]
nc.sids$m.mult <- r$summary.random$idx[1:100, "mean"]
dev.new()
plot(nc.sids$m.uni, nc.sids$m.mult,
main="Posterior mean of the spatial effect", xlab="Uni. post. means"
, ylab="Mult. post. means", xlim=c(-1,1), ylim=c(-1,1))
abline(h=0, col="grey")
abline(v=0, col="grey")
abline(a=0, b=1, col="red")
dev.new()
spplot(nc.sids, c("m.mult", "m.uni"),
main=list(label="Post. mean spatial effect",cex=2))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.