library(sf)
library(tidyverse)
library(spdep)
library(INLA)
library(ggthemes)Spatial and Temporal Analysis Using INLA
Introduction
This document outlines the process of analyzing spatial and tedrivemporal data using Integrated Nested Laplace Approximations (INLA) in R. We aim to explore the effects of several predictors on observed outcomes across different spatial and temporal scales.
Setup
Load Required Libraries
Reading layer `final_sorted' from data source `E:\inla\final_sorted.json' using driver `GeoJSON'
Simple feature collection with 2260 features and 22 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 105.9219 ymin: 19.95672 xmax: 106.5805 ymax: 20.50066
Geodetic CRS: WGS 84theme_Publication <- function(base_size=14) {
library(grid)
library(ggthemes)
(theme_foundation(base_size=base_size)
+ theme(plot.title = element_text(face = "bold",
size = rel(1.2), hjust = 0.5),
text = element_text(),
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour="black", fill=NA, size=0.5),
axis.title = element_text(face = "bold",size = rel(1)),
axis.title.y = element_text(angle=90,vjust =2),
axis.title.x = element_text(vjust = -0.2),
axis.text = element_text(),
axis.line = element_line(colour="black"),
axis.ticks = element_line(),
panel.grid.major = element_line(colour="#f0f0f0"),
panel.grid.minor = element_blank(),
legend.key = element_rect(colour = NA),
plot.margin=unit(c(10,5,5,5),"mm"),
strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
strip.text = element_text(face="bold")
))
}
scale_fill_Publication <- function(...){
library(scales)
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
}
scale_colour_Publication <- function(...){
library(scales)
discrete_scale("colour","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
}Model fitting
Spatial convolution model (1a)
formula1<-adj_Y~1+f(idarea,model="iid")+f(idarea1, model="besag",graph=g)
result1<-inla(formula1,family="poisson",data=data,
E=E,control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result1)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.473, Running = 1.62, Post = 0.31, Total = 2.41
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.142 0.012 0.118 0.142 0.166 0.142 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 50.71 21.15 21.29 46.80 102.98 39.86
Precision for idarea1 17.90 9.70 5.79 15.74 42.72 12.17
Deviance Information Criterion (DIC) ...............: 10968.75
Deviance Information Criterion (DIC, saturated) ....: 2473.09
Effective number of parameters .....................: 164.06
Watanabe-Akaike information criterion (WAIC) ...: 10978.50
Effective number of parameters .................: 162.43
Marginal log-Likelihood: -5738.69
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution model with covariates (model 1b)
formula1b<-adj_Y~1+f(idarea,model="iid")+f(idarea1,model="besag",graph=g) + poor + dens
result1b<-inla(formula1b,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result1b)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.486, Running = 1.49, Post = 0.254, Total = 2.23
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.138 0.020 0.099 0.138 0.177 0.138 0
poor 0.002 0.003 -0.004 0.002 0.008 0.002 0
dens -0.002 0.005 -0.011 -0.002 0.007 -0.002 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 52.23 22.44 21.37 47.98 107.93 40.50
Precision for idarea1 16.96 9.09 5.55 14.95 40.18 11.62
Deviance Information Criterion (DIC) ...............: 10970.14
Deviance Information Criterion (DIC, saturated) ....: 2474.48
Effective number of parameters .....................: 165.46
Watanabe-Akaike information criterion (WAIC) ...: 10980.30
Effective number of parameters .................: 164.00
Marginal log-Likelihood: -5756.43
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + time trend (2)
formula2<-adj_Y~1+f(idarea,model="iid")+f(idarea1,model="besag",graph=g)+idtime
result2<-inla(formula2,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result2)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.273, Running = 1.48, Post = 0.256, Total = 2.01
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.127 0.019 0.089 0.127 0.165 0.127 0
idtime 0.003 0.003 -0.003 0.003 0.008 0.003 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 50.83 21.44 21.02 46.85 103.81 39.80
Precision for idarea1 18.00 9.93 5.74 15.75 43.50 12.08
Deviance Information Criterion (DIC) ...............: 10969.76
Deviance Information Criterion (DIC, saturated) ....: 2474.10
Effective number of parameters .....................: 165.07
Watanabe-Akaike information criterion (WAIC) ...: 10979.56
Effective number of parameters .................: 163.37
Marginal log-Likelihood: -5747.50
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + time trend + cov (2b)
formula2b<-adj_Y~1+f(idarea,model="iid")+f(idarea1,model="besag",graph=g)+idtime + poor + dens
result2b<-inla(formula2b,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result2b)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.261, Running = 1.52, Post = 0.269, Total = 2.05
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.025 0.045 -0.064 0.025 0.114 0.025 0
idtime 0.013 0.005 0.004 0.013 0.022 0.013 0
poor 0.014 0.005 0.004 0.014 0.024 0.014 0
dens -0.001 0.005 -0.010 -0.001 0.008 -0.001 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 49.99 20.76 20.86 46.22 101.05 39.49
Precision for idarea1 18.80 10.70 5.84 16.32 46.42 12.34
Deviance Information Criterion (DIC) ...............: 10966.69
Deviance Information Criterion (DIC, saturated) ....: 2471.03
Effective number of parameters .....................: 166.11
Watanabe-Akaike information criterion (WAIC) ...: 10976.49
Effective number of parameters .................: 164.25
Marginal log-Likelihood: -5761.44
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + temporal random walk 1 (3a)
formula3<-adj_Y~1+f(idarea,model="iid")+f(idarea1, model="besag",graph=g) + f(idtime,model="rw1")
result3<-inla(formula3,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result3)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.346, Running = 1.16, Post = 0.345, Total = 1.85
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.143 0.013 0.118 0.143 0.167 0.143 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
idtime RW1 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 50.76 21.18 21.16 46.88 102.94 39.98
Precision for idarea1 17.99 9.83 5.81 15.77 43.22 12.15
Precision for idtime 26272.40 25863.32 2702.38 18592.22 94889.13 7582.76
Deviance Information Criterion (DIC) ...............: 10969.25
Deviance Information Criterion (DIC, saturated) ....: 2473.59
Effective number of parameters .....................: 165.70
Watanabe-Akaike information criterion (WAIC) ...: 10978.97
Effective number of parameters .................: 163.84
Marginal log-Likelihood: -5740.10
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + temporal random walk 1 + cov (3b)
formula3b<-adj_Y~1+f(idarea,model="iid")+f(idarea1, model="besag",graph=g) + f(idtime,model="rw1") + poor + dens
result3b<-inla(formula3b,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result3b)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.316, Running = 1.26, Post = 0.322, Total = 1.89
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.125 0.024 0.077 0.126 0.170 0.126 0
poor 0.006 0.005 -0.002 0.006 0.016 0.005 0
dens -0.001 0.005 -0.011 -0.001 0.007 -0.001 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
idtime RW1 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 52.33 22.66 21.06 48.07 108.44 40.53
Precision for idarea1 17.16 9.29 5.63 15.07 41.00 11.66
Precision for idtime 19721.36 28299.30 1229.01 11180.65 90792.11 3240.56
Deviance Information Criterion (DIC) ...............: 10968.84
Deviance Information Criterion (DIC, saturated) ....: 2473.18
Effective number of parameters .....................: 167.21
Watanabe-Akaike information criterion (WAIC) ...: 10979.16
Effective number of parameters .................: 165.67
Marginal log-Likelihood: -5756.85
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Bernadinielili
formula4 <- adj_Y ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid") + idtime
result4 <- inla(formula4,
family = "poisson", data = data, E = E,
control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE)
)
summary(result4)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.235, Running = 0.886, Post = 0.323, Total = 1.44
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.128 0.019 0.090 0.128 0.166 0.128 0
idtime 0.003 0.003 -0.003 0.003 0.008 0.003 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 51.32 21.73 21.18 47.28
Precision for idarea (spatial component) 18.56 10.28 5.89 16.23
Precision for idarea1 29529.62 42047.40 2888.09 17206.50
0.975quant mode
Precision for idarea (iid component) 105.03 40.12
Precision for idarea (spatial component) 44.97 12.42
Precision for idarea1 133900.29 7100.70
Deviance Information Criterion (DIC) ...............: 10964.41
Deviance Information Criterion (DIC, saturated) ....: 2468.75
Effective number of parameters .....................: 179.75
Watanabe-Akaike information criterion (WAIC) ...: 10972.23
Effective number of parameters .................: 174.56
Marginal log-Likelihood: -5535.18
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Bernaderlli with covariates (poor and pop dens)
formula4b <- adj_Y ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid") + idtime + poor + dens
result4b <- inla(formula4b,
family = "poisson", data = data, E = E,
control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE)
)
summary(result4b)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.28, Running = 1.21, Post = 0.531, Total = 2.02
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.026 0.046 -0.064 0.026 0.116 0.026 0
idtime 0.013 0.005 0.003 0.013 0.022 0.013 0
poor 0.014 0.005 0.004 0.014 0.024 0.014 0
dens -0.001 0.005 -0.010 -0.001 0.008 -0.001 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 49.97 21.20 20.29 46.11
Precision for idarea (spatial component) 20.23 12.36 5.86 17.23
Precision for idarea1 30246.99 43666.78 2871.27 17466.13
0.975quant mode
Precision for idarea (iid component) 102.14 39.18
Precision for idarea (spatial component) 52.45 12.60
Precision for idarea1 138261.90 7085.42
Deviance Information Criterion (DIC) ...............: 10961.46
Deviance Information Criterion (DIC, saturated) ....: 2465.80
Effective number of parameters .....................: 180.41
Watanabe-Akaike information criterion (WAIC) ...: 10969.38
Effective number of parameters .................: 175.17
Marginal log-Likelihood: -5549.08
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + temporal random walk + idd space-time interaction
formula5<-adj_Y~1+f(idarea,model="iid")+f(idarea1,model="besag",graph=g)+ f(idtime,model="rw1")+f(id2,model="iid")
result5<-inla(formula5,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result5)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.294, Running = 2.38, Post = 0.441, Total = 3.12
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.142 0.012 0.118 0.142 0.167 0.142 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
idtime RW1 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 50.96 21.59 20.94 46.97 104.28 39.86
Precision for idarea1 17.98 9.95 5.74 15.72 43.56 12.05
Precision for idtime 26279.80 25730.12 2784.30 18665.90 94527.88 7794.15
Precision for id2 21165.24 23739.03 1260.21 13625.87 84157.26 3358.56
Deviance Information Criterion (DIC) ...............: 10968.38
Deviance Information Criterion (DIC, saturated) ....: 2472.72
Effective number of parameters .....................: 167.36
Watanabe-Akaike information criterion (WAIC) ...: 10978.54
Effective number of parameters .................: 165.68
Marginal log-Likelihood: -5740.12
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Spatial convolution + temporal random walk + idd space-time interaction + cov
formula5b<-adj_Y~1+f(idarea,model="iid")+f(idarea1,model="besag",graph=g)+ f(idtime,model="rw1")+f(id2,model="iid") + poor + dens
result5b<-inla(formula5b,family="poisson",data=data,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result5b)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.287, Running = 2.9, Post = 0.622, Total = 3.81
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.125 0.024 0.077 0.125 0.170 0.126 0
poor 0.006 0.005 -0.002 0.006 0.016 0.005 0
dens -0.001 0.005 -0.010 -0.001 0.007 -0.001 0
Random effects:
Name Model
idarea IID model
idarea1 Besags ICAR model
idtime RW1 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for idarea 51.91 22.50 20.94 47.66 107.71 40.15
Precision for idarea1 17.54 9.75 5.55 15.32 42.60 11.71
Precision for idtime 21350.09 35124.02 984.50 10940.82 106206.96 2491.70
Precision for id2 21116.11 23227.88 1296.99 13759.43 82768.12 3482.56
Deviance Information Criterion (DIC) ...............: 10968.44
Deviance Information Criterion (DIC, saturated) ....: 2472.78
Effective number of parameters .....................: 169.41
Watanabe-Akaike information criterion (WAIC) ...: 10978.76
Effective number of parameters .................: 167.59
Marginal log-Likelihood: -5756.80
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Analyzing M4b - best fitted
result4b$summary.fixed mean sd 0.025quant 0.5quant 0.975quant
(Intercept) 0.0257043857 0.045998475 -0.064456089 0.0256889120 0.115953101
idtime 0.0128028233 0.004777249 0.003427921 0.0128051095 0.022164647
poor 0.0137405427 0.005141105 0.003655433 0.0137418013 0.023818481
dens -0.0007320154 0.004517961 -0.009763303 -0.0006774755 0.007987909
mode kld
(Intercept) 0.0256888776 4.478116e-11
idtime 0.0128051023 7.456053e-11
poor 0.0137418058 4.872449e-11
dens -0.0006778031 3.850365e-08result4b$summary.hyperpar mean sd 0.025quant
Precision for idarea (iid component) 49.97340 21.20415 20.291246
Precision for idarea (spatial component) 20.23007 12.35707 5.860052
Precision for idarea1 30246.98459 43666.78032 2871.265844
0.5quant 0.975quant mode
Precision for idarea (iid component) 46.10784 102.13994 39.17711
Precision for idarea (spatial component) 17.22799 52.44853 12.59474
Precision for idarea1 17466.13065 138261.89763 7085.41877Sensitivity analysis
formula4b_p3 <- adj_Y ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid", param=c(0.1, 0.5)) + idtime + poor + dens
result4b_p3 <- inla(formula4b_p3,
family = "poisson", data = data, E = E,
control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE)
)
summary(result4b_p3)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.298, Running = 0.937, Post = 0.498, Total = 1.73
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.057 0.055 -0.051 0.057 0.165 0.057 0
idtime 0.009 0.008 -0.006 0.009 0.024 0.009 0
poor 0.012 0.006 0.000 0.012 0.024 0.012 0
dens -0.006 0.005 -0.016 -0.006 0.005 -0.006 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 17.16 3.79 10.88 16.77
Precision for idarea (spatial component) 62.34 65.72 9.52 42.99
Precision for idarea1 158.09 15.61 129.28 157.43
0.975quant mode
Precision for idarea (iid component) 25.73 16.02
Precision for idarea (spatial component) 235.12 22.61
Precision for idarea1 190.67 156.33
Deviance Information Criterion (DIC) ...............: 10988.31
Deviance Information Criterion (DIC, saturated) ....: 2492.65
Effective number of parameters .....................: 322.78
Watanabe-Akaike information criterion (WAIC) ...: 10980.14
Effective number of parameters .................: 278.05
Marginal log-Likelihood: -5746.85
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')formula4b_bym_iid <- adj_Y ~ f(idarea, model = "bym", graph = g, hyper = list( prec.spatial = list(prior = "loggamma", param = c(0.01, 0.01)), prec.unstruct = list(prior = "loggamma", param = c(0.01, 0.01)) )) + f(idarea1, idtime, model = "iid", hyper = list( prec = list(prior = "loggamma", param = c(0.01, 0.01)) )) + idtime + poor + dens
result4b_bym_idd <- inla(formula4b_bym_iid, family = "poisson", data = data, E = E, control.predictor = list(compute = TRUE), control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE) )
summary(result4b_bym_idd)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.293, Running = 1.01, Post = 0.502, Total = 1.8
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.034 0.048 -0.061 0.034 0.128 0.034 0
idtime 0.012 0.005 0.001 0.012 0.022 0.012 0
poor 0.013 0.005 0.003 0.013 0.024 0.013 0
dens -0.002 0.005 -0.011 -0.002 0.007 -0.002 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 40.24 13.54 19.66 38.24
Precision for idarea (spatial component) 23.20 13.48 7.29 19.98
Precision for idarea1 1585.59 266.25 1126.61 1563.53
0.975quant mode
Precision for idarea (iid component) 72.33 34.57
Precision for idarea (spatial component) 58.30 15.00
Precision for idarea1 2171.82 1520.06
Deviance Information Criterion (DIC) ...............: 10955.38
Deviance Information Criterion (DIC, saturated) ....: 2459.72
Effective number of parameters .....................: 229.83
Watanabe-Akaike information criterion (WAIC) ...: 10957.83
Effective number of parameters .................: 212.10
Marginal log-Likelihood: -5580.05
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')formula4b_l <- adj_Y ~ f(idarea, model = "bym",hyper = list( prec.spatial = list(prior = "loggamma", param = c(2, 0.05)), prec.unstruct = list(prior = "loggamma", param = c(2, 0.5)) ), graph = g) + f(idarea1, idtime, model = "iid", hyper = list( prec = list(prior = "loggamma", param = c(2, 0.5)) )) + idtime + poor + dens
result4b_l <- inla(formula4b_l, family = "poisson", data = data, E = E, control.predictor = list(compute = TRUE), control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE))
summary(result4b_l)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.512, Running = 0.97, Post = 0.432, Total = 1.91
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.058 0.056 -0.052 0.058 0.168 0.058 0
idtime 0.009 0.008 -0.006 0.009 0.024 0.009 0
poor 0.012 0.006 0.000 0.012 0.024 0.012 0
dens -0.006 0.006 -0.017 -0.006 0.005 -0.006 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 14.15 2.36 10.06 13.96
Precision for idarea (spatial component) 40.67 22.29 12.89 35.70
Precision for idarea1 158.68 15.57 130.02 157.99
0.975quant mode
Precision for idarea (iid component) 19.33 13.59
Precision for idarea (spatial component) 97.72 27.44
Precision for idarea1 191.25 156.79
Deviance Information Criterion (DIC) ...............: 10984.04
Deviance Information Criterion (DIC, saturated) ....: 2488.38
Effective number of parameters .....................: 329.71
Watanabe-Akaike information criterion (WAIC) ...: 10973.90
Effective number of parameters .................: 281.88
Marginal log-Likelihood: -5732.18
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')formula4b_l3 <- adj_Y ~ f(idarea, model = "bym",hyper = list( prec.spatial = list(prior = "loggamma", param = c(1, 0.0001)), prec.unstruct = list(prior = "loggamma", param = c(1, 0.0001)) ), graph = g) + f(idarea1, idtime, model = "iid", hyper = list( prec = list(prior = "loggamma", param = c(1, 0.0001)) )) + idtime + poor + dens
result4b_l3 <- inla(formula4b_l3, family = "poisson", data = data, E = E, control.predictor = list(compute = TRUE), control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE))
summary(result4b_l3)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " keep = keep, working.directory = working.directory,
silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
debug, .parent.frame = .parent.frame)" )
Time used:
Pre = 0.267, Running = 1.03, Post = 0.43, Total = 1.73
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.026 0.046 -0.065 0.026 0.117 0.026 0
idtime 0.013 0.005 0.003 0.013 0.022 0.013 0
poor 0.014 0.005 0.004 0.014 0.024 0.014 0
dens -0.001 0.005 -0.010 -0.001 0.008 -0.001 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 50.98 22.35 20.29 46.74
Precision for idarea (spatial component) 20.27 12.63 5.61 17.20
Precision for idarea1 12961.61 10389.48 2665.86 10096.68
0.975quant mode
Precision for idarea (iid component) 106.43 39.25
Precision for idarea (spatial component) 53.16 12.40
Precision for idarea1 40559.45 6268.64
Deviance Information Criterion (DIC) ...............: 10960.06
Deviance Information Criterion (DIC, saturated) ....: 2464.40
Effective number of parameters .....................: 183.51
Watanabe-Akaike information criterion (WAIC) ...: 10968.08
Effective number of parameters .................: 177.92
Marginal log-Likelihood: -5552.46
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')id time effect:
ggplot(as.data.frame(result4b_bym_idd$marginals.fixed$idtime)) +
geom_line(aes(x=x,y=y))
plot_combined_marginal_effects_with_quantiles <- function(result_list, effect_name, model_names) {
# Create an empty data frame to store combined data
combined_data <- data.frame(x = numeric(), y = numeric(), Model = factor(), lquant = numeric(), uquant = numeric())
# Loop through the list of model results
for (i in seq_along(result_list)) {
# Extract the marginal effects for the specified variable
effect <- result_list[[i]]$marginals.fixed[[effect_name]]
# Calculate the quantiles for the shading
lquant <- inla.qmarginal(0.025, effect)
uquant <- inla.qmarginal(0.975, effect)
# Create a data frame from the smoothed marginal
df_effect <- data.frame(inla.smarginal(effect))
df_effect$Model <- model_names[i]
df_effect$lquant <- lquant
df_effect$uquant <- uquant
# Combine with the main data frame
combined_data <- rbind(combined_data, df_effect)
}
# Define line types and widths for better differentiation
line_types <- c("solid", "dashed", "dotdash", "longdash")
line_widths <- c(1.2, 1.2, 1.2, 1.2)
# Generate the plot
p <- ggplot(combined_data, aes(x, y, color = Model)) +
geom_line(aes(linetype = Model), size = 1.2) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_colour_Publication() +
scale_linetype_manual(values = line_types) +
theme_Publication(base_size = 14) +
theme(legend.position = "bottom", legend.title = element_blank()) +
labs(title = NULL,
x = NULL, y = "Density")
return(p)
}
# Usage
# Assuming you have result objects like result4b_bym_idd, result4b_l3, result4b_l, result4b
result_list <- list(result4b, result4b_bym_idd, result4b_l)
model_names <- c("Default priors, log-gamma (1, 0.0005)","Best-fitted, log-gamma (0.01, 0.01)", "log-gamma (2, 0.5)")
plot_combined_marginal_effects_with_quantiles(result_list, "poor", model_names)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Attaching package: 'scales'The following object is masked from 'package:purrr':
discardThe following object is masked from 'package:readr':
col_factorWarning: The `scale_name` argument of `discrete_scale()` is deprecated as of ggplot2
3.5.0.Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
x <- 1:10 # This represents the years in a relative sense.
years <- 2013:2022 # Actual years corresponding to each x value
# Assuming model.par$summary.fixed[2,1] and others are your model parameters
beta_main <- result4b_bym_idd$summary.fixed[2,1]
beta_lower <- result4b_bym_idd$summary.fixed[2,3]
beta_upper <- result4b_bym_idd$summary.fixed[2,5]
# Create a data frame for plotting
df <- data.frame(
Year = years,
Main = beta_main * x,
Lower = beta_lower * x,
Upper = beta_upper * x
)
# Melting data frame for easier plotting with ggplot2
df_long <- tidyr::pivot_longer(df, cols = c("Main", "Lower", "Upper"), names_to = "LineType", values_to = "Value")
time_eff <- ggplot(df_long, aes(x = Year, y = Value, group = LineType, color = LineType)) +
geom_line(aes(linetype = LineType)) + # Add lines with different types
scale_linetype_manual(values = c("Main" = "solid", "Lower" = "dashed", "Upper" = "dashed")) + # Customize line types
scale_color_manual(values = c("Main" = "black", "Lower" = "red", "Upper" = "blue")) + # Customize line colors if needed
labs(x = "", y = "", title = "") + # Add labels
theme_Publication() + # Apply the publication theme
theme(
legend.title = element_blank(), # Remove legend title
legend.position = "none" # Remove legend if not needed
) + scale_x_continuous(breaks = 2013:2022, labels = as.character(2013:2022))
time_eff
plot_combined_marginal_effects_with_quantiles <- function(result_list, effect_name, model_names) {
# Create an empty data frame to store combined data
combined_data <- data.frame(x = numeric(), y = numeric(), Model = factor(), lquant = numeric(), uquant = numeric())
# Loop through the list of model results
for (i in seq_along(result_list)) {
# Extract the marginal effects for the specified variable
effect <- result_list[[i]]$marginals.fixed[[effect_name]]
# Calculate the quantiles for the shading
lquant <- inla.qmarginal(0.025, effect)
uquant <- inla.qmarginal(0.975, effect)
# Create a data frame from the smoothed marginal
df_effect <- data.frame(inla.smarginal(effect))
df_effect$Model <- model_names[i]
df_effect$lquant <- lquant
df_effect$uquant <- uquant
# Combine with the main data frame
combined_data <- rbind(combined_data, df_effect)
}
# Define line types and widths for better differentiation
line_types <- c("solid", "dashed", "dotdash", "longdash")
line_widths <- c(1.2, 1.2, 1.2, 1.2)
# Generate the plot
p <- ggplot(combined_data, aes(x, y, color = Model)) +
geom_line(aes(linetype = Model), size = 1.2) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_colour_Publication() +
scale_linetype_manual(values = line_types) +
theme_Publication(base_size = 14) +
theme(legend.position = "bottom", legend.title = element_blank()) +
labs(title = NULL,
x = NULL, y = "Probability density")
return(p)
}
# Usage
# Assuming you have result objects like result4b_bym_idd, result4b_l3, result4b_l, result4b
result_list <- list(result4b_bym_idd, result4b, result4b_l)
model_names <- c("Best-fitted, log-gamma (0.01, 0.01)", "Default priors, log-gamma (1, 0.0005)", "log-gamma (2, 0.5)")
plot_combined_marginal_effects_with_quantiles(result_list, "poor", model_names)