# NOT RUN {
# }
# NOT RUN {
if(require(INLA, quietly = TRUE)) {
require(sp)
require(spdep)
data(CV)
W <- as(nb2mat(CV.nb, style = "B"), "Matrix")
#Data (two diseases only)
d <- list(OBS = c(CV$Obs.Cirrhosis, CV$Obs.Lung),
EXP = c(CV$Exp.Cirrhosis, CV$Exp.Lung))
# Index for latent effect
d$idx <- 1:length(d$OBS)
k <- 2 #Number of diseases
# Linear constraint for models
A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
e = rep(0, k)
# Two independent ICAR models
#model <- inla.rgeneric.define(inla.rgeneric.indep.IMCAR.model,
# k = k, W = W)
model <- inla.INDIMCAR.model(k = k, W = W)
r.simcar <- try(
inla(OBS ~ 1 + f(idx, model = model, extraconstr = list(A = as.matrix(A), e = e)),
data = d, E = EXP, family = "poisson",
# To run faster, REMOVE in real applications
control.mode = list(theta = c(1.4, 2.1), restart = TRUE),
control.predictor = list(compute = TRUE))
)
summary(r.simcar)
# IMCAR model
#model <- inla.rgeneric.define(inla.rgeneric.IMCAR.model,
# k = k, W = W, alpha.min = 0, alpha.max = 1)
model <- inla.IMCAR.model(k = k, W = W)
r.imcar <- try(
inla(OBS ~ 1 + f(idx, model = model, extraconstr = list(A = as.matrix(A), e = e)),
data = d, E = EXP, family = "poisson",
# To run faster, REMOVE in real applications
control.mode = list(theta = c(1.77, 2.01, 0.93),
restart = TRUE),
control.compute = list(config = TRUE),
control.predictor = list(compute = TRUE))
)
summary(r.imcar)
# Transform parameters
summary.post <- inla.MCAR.transform(r.imcar, k = k)
# Posterior of variance matrix
summary.post$VAR.p # Using point estimates
summary.post$VAR.m # Using posterior sampling
} #if(require(INLA))
# }
Run the code above in your browser using DataLab