# NOT RUN {
# }
# NOT RUN {
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))
}
# }
# NOT RUN {
# }
Run the code above in your browser using DataCamp Workspace