indepimcar: MCAR(1, Lambda): Intrinsic multivariate CAR latent effect...

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

Description

Multivariate generalization of the intrinsic conditional autorregresive model. No correlation parameters are considered between the different diseases, so the matrix which models the variability between diseases will be a diagonal matrix.

Usage

1
2
3

Arguments

...

Arguments to be passed to 'inla.rgeneric.define'.

cmd

Arguments used by latent effects defined using the 'rgeneric' latent effect.

theta

Vector of hyperparameters.

Details

This function is used to define a latent effect that is a multivariate spatial effect with a intrinsic conditional autorregresive distribution and a diagonal 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:

This model is defined using the 'f()' function and an index in order to identify the spatial areas. See the example.

Value

This is used internally by the 'INLA::inla()'.

Prior distributions of the hyperparameters

The hyperparamenters of this lattent effect are the marginal precisions of each disease. So the total number of hyperpameters is equal to the number of diseases.

References

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

Examples

  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
if (require("INLA", quietly = TRUE)) {
require(spdep)
require(spData)
require(rgdal)

## Independent IMCAR model with 2 diseases

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

# Model parameters: k and W
k <- 2 * n.rep #Number of diseases

#Define independent IMCAR model
model <- inla.rgeneric.define(inla.rgeneric.indep.IMCAR.model, debug = FALSE,
                             k = k,
                             W = W)

# Matrices for sum-to-zero constraints
A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
e  = rep(0, k)

#Fit multivariate 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.predictor = list(compute = TRUE))

summary(r)

# Transformed parameters
r.hyperpar <- inla.MCAR.transform(r, k = 2, model = "INDIMCAR")
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, alpha

marg.tau1 <- inla.tmarginal(
 function(x) exp(x),
 r$marginals.hyperpar[[1]])

marg.tau2 <- inla.tmarginal(
 function(x) exp(x),
 r$marginals.hyperpar[[2]])

oldpar <- par(mfrow = c(2, 1))
plot(marg.tau1, main = "tau1", type = "l")
plot(marg.tau2, main = "tau2", type = "l")

par(oldpar)

## Running UNIVARIATE MODEL

#Real data
d.uni <- list(OBS = nc.sids$SID74,
             NWPROP = nc.sids$NWPROP74,
             EXP = nc.sids$EXP74)
d.uni$idx <- 1:length(d.uni$OBS)

#Fit model
r.uni <- inla(OBS ~ 1 + f(idx, model = "besag", graph = W),
             data = d.uni, 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.

spplot(nc.sids, c("FITTED74", "FITTED74.uni"),
      main=list(label="Relative risk estimation",cex=2))

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"]

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

spplot(nc.sids, c("m.mult", "m.uni"),
      main=list(label="Post. mean spatial effect",cex=2))

}

INLAMSM documentation built on June 4, 2021, 9:07 a.m.