# packages
library(INLA)
#> Loading required package: Matrix
#> Loading required package: sp
#> Loading required package: parallel
#> Loading required package: foreach
#> This is INLA_20.02.19 built 2020-02-19 06:22:22 UTC.
#> See www.r-inla.org/contact-us for how to get help.
library(igraph)
#>
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#>
#> decompose, spectrum
#> The following object is masked from 'package:base':
#>
#> union
# simulate data
set.seed(02032020)
n <- 100
a <- 0.5
b <- 1.5
x1 <- rnorm(n, sd = 0.5)
eta.z <- -a - b*x1
z <- rbinom(n, 1, inla.link.logit(eta.z, inverse=TRUE))
n.y <- sum(z)
x2 <- rnorm(n.y, sd = 0.5)
eta.y <- a + b*x2
lambda <- exp(eta.y)
y <- rpois(n.y, lambda)
is.zero <- (y == 0)
while(sum(is.zero) > 0) {
y[is.zero] = rpois(sum(is.zero), lambda[is.zero])
is.zero = (y == 0)
}
# Create the appropriate Y matrix
Y <- matrix(NA, n + n.y, 2)
Y[1:n, 1] <- z
Y[n + 1:n.y, 2] <- y
form = Y ~ 0 + mu.z + mu.y + cov.z + cov.y
ldat <- list(
Y = Y,
mu.z = rep(1:0, c(n, n.y)),
mu.y = rep(0:1, c(n, n.y)),
cov.z = c(x1, rep(NA, n.y)),
cov.y = c(rep(NA, n), x2)
)
# Estimate the model
summary(inla(
form,
data = ldat,
family=c('binomial', 'zeroinflatedpoisson0'),
control.family = list(
list(),
list(
hyper = list(
prob = list(
initial = -20,
fixed = TRUE
)
)
)
)
))
#>
#> Call:
#> c("inla(formula = form, family = c(\"binomial\",
#> \"zeroinflatedpoisson0\"), ", " data = ldat, control.family =
#> list(list(), list(hyper = list(prob = list(initial = -20, ", " fixed =
#> TRUE)))))")
#> Time used:
#> Pre = 1.31, Running = 0.232, Post = 0.0339, Total = 1.58
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> mu.z -0.134 0.214 -0.556 -0.133 0.286 -0.132 0
#> mu.y 0.586 0.135 0.309 0.591 0.837 0.600 0
#> cov.z -1.421 0.427 -2.304 -1.406 -0.629 -1.374 0
#> cov.y 1.219 0.208 0.810 1.219 1.628 1.219 0
#>
#> Expected number of effective parameters(stdev): 4.00(0.00)
#> Number of equivalent replicates : 36.50
#>
#> Marginal log-Likelihood: -147.31
# Define the spatial structure
random_graph <- sample_k_regular(n, 3)
par(mar = rep(0, 4))
plot(random_graph)
spatial_structure <- as_adjacency_matrix(random_graph)
# Add the spatial structure to the model
ids_z <- seq(1, n)
ids_y <- ids_z[z != 0]
# Redefine formula and data
form = Y ~ 0 + mu.z + mu.y + cov.z + cov.y +
f(spatial_z, model = "bym", graph = spatial_structure) +
f(spatial_y, model = "bym", graph = spatial_structure)
ldat <- list(
Y = Y,
mu.z = rep(1:0, c(n, n.y)),
mu.y = rep(0:1, c(n, n.y)),
cov.z = c(x1, rep(NA, n.y)),
cov.y = c(rep(NA, n), x2),
spatial_z = c(ids_z, rep(NA, n.y)),
spatial_y = c(rep(NA, n), ids_y)
)
summary(inla(
form,
data = ldat,
family=c('binomial', 'zeroinflatedpoisson0'),
control.family = list(
list(),
list(
hyper = list(
prob = list(
initial = -20,
fixed = TRUE
)
)
)
)
))
#>
#> Call:
#> c("inla(formula = form, family = c(\"binomial\",
#> \"zeroinflatedpoisson0\"), ", " data = ldat, control.family =
#> list(list(), list(hyper = list(prob = list(initial = -20, ", " fixed =
#> TRUE)))))")
#> Time used:
#> Pre = 1.76, Running = 0.831, Post = 0.0699, Total = 2.66
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> mu.z -0.134 0.214 -0.556 -0.133 0.286 -0.132 0
#> mu.y 0.585 0.135 0.307 0.589 0.836 0.599 0
#> cov.z -1.422 0.427 -2.305 -1.406 -0.629 -1.375 0
#> cov.y 1.221 0.209 0.810 1.220 1.632 1.220 0
#>
#> Random effects:
#> Name Model
#> spatial_z BYM model
#> spatial_y BYM model
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant
#> Precision for spatial_z (iid component) 2580.00 2788.60 244.56 1752.12
#> Precision for spatial_z (spatial component) 2674.66 2924.62 260.32 1804.84
#> Precision for spatial_y (iid component) 2229.05 2247.47 198.56 1570.86
#> Precision for spatial_y (spatial component) 2214.96 2254.06 188.98 1552.15
#> 0.975quant mode
#> Precision for spatial_z (iid component) 9920.81 670.62
#> Precision for spatial_z (spatial component) 10332.83 714.34
#> Precision for spatial_y (iid component) 8178.17 555.86
#> Precision for spatial_y (spatial component) 8180.32 530.37
#>
#> Expected number of effective parameters(stdev): 4.17(0.116)
#> Number of equivalent replicates : 35.03
#>
#> Marginal log-Likelihood: -49.82Created on 2020-03-02 by the reprex package (v0.3.0)