Mmodel: M-model: Proper multivariate CAR latent effect with a...
In INLAMSM: Multivariate Spatial Models with 'INLA'

Description Usage Arguments Details Value Prior distributions of the hyperparameters References Examples

Multivariate generalization of the proper conditional autorregresive model with one common correlation parameter. This model is performed using the M-model aproximation of Rocamora et. al. (2015).

1
2
3

inla.rgeneric.Mmodel.model(cmd, theta)

inla.Mmodel.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 based on the M-model aproximation of Rocamora et. al. (2015) in which θ is modelled as a product of a Φ \cdot M where the colums of Φ are modeled independently with a proper conditional autorregresive distribution with a different spatial autocorrelation parameter for each disease and M is a square matrix which introduce de dependence between the diseases. 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.
alpha.min Minimum value of the spatial autocorrelation parameter.
alpha.max Maximum value of the spatial autocorrelation parameter.

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 common spatial autocorrelation parameters (one for each disease) and the entries of the M matrix (considered all as a random effects).

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).

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

# 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)

#Parameters of the Mmodel
k <- 2
alpha.min <- 0
alpha.max <- 1



model <- inla.rgeneric.define(inla.rgeneric.Mmodel.model, debug = FALSE,
                              k = k, W = W, alpha.min = alpha.min,
                              alpha.max = alpha.max)

r.Mmodel <- inla(OBS ~ -1 + f(idx, model = model), data = d, E = EXP,
  family = "poisson", control.predictor = list(compute = TRUE))

nc.sids$Model1 <- r.Mmodel$summary.random$idx[1:100, "mean"]
nc.sids$Model2 <- r.Mmodel$summary.random$idx[100 + 1:100, "mean"]

spplot(nc.sids, c("Model1", "Model2"))

nc.sids$Fit1 <- r.Mmodel$summary.fitted[1:100, "mean"]
nc.sids$Fit2 <- r.Mmodel$summary.fitted[100 + 1:100, "mean"]

spplot(nc.sids, c("Fit1", "SMR74", "Fit2", "SMR79"))


## 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", "Fit1", "FITTED74.uni"))
spplot(nc.sids, c("Fit1", "FITTED74.uni"),
       main=list(label="Relative risk estimation",cex=2))
dev.new()
plot(nc.sids$FITTED74.uni, nc.sids$Fit1, main="Relative Risk estimations",
     xlab="Univariate RR estimations"
     , ylab="Multivariate RR estimations")#, xlim=c(0.5, 2.5), ylim=c(0, 2))
abline(h=0, col="grey")
abline(v=0, col="grey")
abline(a=0, b=1, col="red")

}