library(sf)
library(tidyverse)
library(spdep)
library(INLA)
library(ggthemes)
library(knitr)Spatial and Temporal Analysis of TB in Nam Dinh Using INLA
Introduction
This document outlines the process of analyzing spatial and temporal 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 `merged_df1' from data source `E:\inla\merged_df1.geojson' using driver `GeoJSON'
Simple feature collection with 2260 features and 24 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=12) {
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, linewidth=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")), ...)
}Calculate TB notifications by time
notif_data <- gdf %>%
group_by(year) %>%
summarise(
Total_Observed = sum(adjusted_observed),
Total_Population = sum(pop)
) %>%
mutate(notif = (Total_Observed / Total_Population) * 1e5)
notif_dataSimple feature collection with 10 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 105.9219 ymin: 19.95672 xmax: 106.5805 ymax: 20.50066
Geodetic CRS: WGS 84
# A tibble: 10 × 5
year Total_Observed Total_Population geometry notif
* <int> <dbl> <int> <POLYGON [°]> <dbl>
1 2013 1488 1846126 ((106.1746 19.99359, 106.1749 19… 80.6
2 2014 1457 1851759 ((106.1746 19.99359, 106.1749 19… 78.7
3 2015 1535 1796891 ((106.1746 19.99359, 106.1749 19… 85.4
4 2016 1556 1792452 ((106.1746 19.99359, 106.1749 19… 86.8
5 2017 1753 1788225 ((106.1746 19.99359, 106.1749 19… 98.0
6 2018 1654 1784206 ((106.1746 19.99359, 106.1749 19… 92.7
7 2019 1688 1780865 ((106.1746 19.99359, 106.1749 19… 94.8
8 2020 1369 1780333 ((106.1746 19.99359, 106.1749 19… 76.9
9 2021 1111 1836268 ((106.1746 19.99359, 106.1749 19… 60.5
10 2022 1276 1830794 ((106.1746 19.99359, 106.1749 19… 69.7ggplot(notif_data, aes(x = year, y = notif)) +
geom_line() +
geom_point() +
theme_Publication() +
labs(
title = "",
x = "Year",
y = "TB case notifications per 100,000"
)
Moran I
weights <- poly2nb(subset_gdf, queen = TRUE)
weights <- nb2listw(weights, style = "W", zero.policy = TRUE)
gdf$adjusted_observed[gdf$adjusted_observed==0] <- NA
gdf$SMR <- gdf$adjusted_observed / gdf$Eunique_years <- unique(gdf$year)
# Initialize an empty list to store results
moran_results <- list()
# Loop through each year and calculate Moran's I
for (current_year in unique_years) {
# Subset data for the specific year
year_data <- dplyr::filter(gdf, year == current_year)
# Check if year_data is empty after subsetting
if (nrow(year_data) == 0) {
cat(sprintf("Moran's I for year %d: No data for this year\n", current_year))
moran_results[[as.character(current_year)]] <- list(year = current_year, moran = NA, p_value = NA, warning = "No data for this year")
next
}
# Convert 'SMR' to numeric
year_data$SMR <- as.numeric(as.character(year_data$SMR))
# Replace NA values in SMR with 0
year_data$SMR[is.na(year_data$SMR)] <- 0
# Ensure the order of SMR values matches the spatial units in the weights matrix
if (length(year_data$SMR) != length(weights$neighbours)) {
cat(sprintf("Moran's I for year %d: Mismatch in data length and weights matrix\n", current_year))
moran_results[[as.character(current_year)]] <- list(year = current_year, moran = NA, p_value = NA, warning = "Mismatch in data length and weights matrix")
next
}
# Perform Moran's I test
moran_test <- moran.test(year_data$SMR, weights, zero.policy = TRUE)
cat(sprintf("Moran's I for year %d is %f with a normal p-value of %f\n", current_year, moran_test$estimate, moran_test$p.value))
# Store results in the list
moran_results[[as.character(current_year)]] <- list(year = current_year, moran = moran_test$estimate, p_value = moran_test$p.value)
}Moran's I for year 2013 is 0.094222 with a normal p-value of 0.007713
Moran's I for year 2013 is -0.004444 with a normal p-value of 0.007713
Moran's I for year 2013 is 0.001659 with a normal p-value of 0.007713
Moran's I for year 2014 is 0.117018 with a normal p-value of 0.001499
Moran's I for year 2014 is -0.004444 with a normal p-value of 0.001499
Moran's I for year 2014 is 0.001675 with a normal p-value of 0.001499
Moran's I for year 2015 is 0.081503 with a normal p-value of 0.017892
Moran's I for year 2015 is -0.004444 with a normal p-value of 0.017892
Moran's I for year 2015 is 0.001676 with a normal p-value of 0.017892
Moran's I for year 2016 is 0.090854 with a normal p-value of 0.010018
Moran's I for year 2016 is -0.004444 with a normal p-value of 0.010018
Moran's I for year 2016 is 0.001679 with a normal p-value of 0.010018
Moran's I for year 2017 is 0.070746 with a normal p-value of 0.033001
Moran's I for year 2017 is -0.004444 with a normal p-value of 0.033001
Moran's I for year 2017 is 0.001673 with a normal p-value of 0.033001
Moran's I for year 2018 is 0.059970 with a normal p-value of 0.057819
Moran's I for year 2018 is -0.004444 with a normal p-value of 0.057819
Moran's I for year 2018 is 0.001676 with a normal p-value of 0.057819
Moran's I for year 2019 is 0.129636 with a normal p-value of 0.000533
Moran's I for year 2019 is -0.004444 with a normal p-value of 0.000533
Moran's I for year 2019 is 0.001678 with a normal p-value of 0.000533
Moran's I for year 2020 is 0.069116 with a normal p-value of 0.035978
Moran's I for year 2020 is -0.004444 with a normal p-value of 0.035978
Moran's I for year 2020 is 0.001671 with a normal p-value of 0.035978
Moran's I for year 2021 is 0.043341 with a normal p-value of 0.122065
Moran's I for year 2021 is -0.004444 with a normal p-value of 0.122065
Moran's I for year 2021 is 0.001683 with a normal p-value of 0.122065
Moran's I for year 2022 is 0.051127 with a normal p-value of 0.086796
Moran's I for year 2022 is -0.004444 with a normal p-value of 0.086796
Moran's I for year 2022 is 0.001668 with a normal p-value of 0.086796# Assuming moran_results is a list of lists with each year's result
# Convert the list to a data frame
results_df <- do.call(rbind, lapply(names(moran_results), function(year) {
moran <- moran_results[[year]]$moran
data.frame(
Year = as.integer(year),
Moran_I = moran[1], # Moran I statistic
Expectation = moran[2], # Expected value
Variance = moran[3], # Variance of Moran's I
P_Value = moran_results[[year]]$p_value
)
}))
### 2. Presenting the Data
# Use kable from knitr to create a nice table
knitr::kable(results_df, format = "html", table.attr = "style='width:100%;'",
col.names = c("Year", "Moran I", "Expectation", "Variance", "P-Value"))| Year | Moran I | Expectation | Variance | P-Value | |
|---|---|---|---|---|---|
| Moran I statistic | 2013 | 0.0942219 | -0.0044444 | 0.0016593 | 0.0077132 |
| Moran I statistic1 | 2014 | 0.1170183 | -0.0044444 | 0.0016749 | 0.0014993 |
| Moran I statistic2 | 2015 | 0.0815028 | -0.0044444 | 0.0016760 | 0.0178917 |
| Moran I statistic3 | 2016 | 0.0908540 | -0.0044444 | 0.0016791 | 0.0100176 |
| Moran I statistic4 | 2017 | 0.0707463 | -0.0044444 | 0.0016728 | 0.0330005 |
| Moran I statistic5 | 2018 | 0.0599697 | -0.0044444 | 0.0016761 | 0.0578186 |
| Moran I statistic6 | 2019 | 0.1296363 | -0.0044444 | 0.0016785 | 0.0005326 |
| Moran I statistic7 | 2020 | 0.0691162 | -0.0044444 | 0.0016712 | 0.0359782 |
| Moran I statistic8 | 2021 | 0.0433409 | -0.0044444 | 0.0016832 | 0.1220651 |
| Moran I statistic9 | 2022 | 0.0511273 | -0.0044444 | 0.0016678 | 0.0867960 |
gdf$year <- as.integer(gdf$year)
data1 <- gdf %>%
mutate(SMR = ifelse(is.na(SMR), 0, SMR)) %>% # Replace NA values with 0
group_by(year) %>%
do({
.data <- .
local_moran_results <- localmoran(.data$SMR, weights)
# Transform results into a dataframe to rejoin with original data
data.frame(.data,
local_I = local_moran_results[, "Ii"],
p_value = 2 * pnorm(-abs(local_moran_results[, "Z.Ii"]), lower.tail = TRUE), # two-sided p-value
hotspot_type = attr(local_moran_results, "quadr")[["pysal"]] # assuming 'pysal' gives the most consistent results
)
}) %>%
ungroup()# Visualization using ggplot2 with facets for each year
merged_data <- gdf %>%
# Step 1: Remove the 'observed' column from gdf before merging
select(-observed, -hotspot_type) %>%
# Step 2: Merge with selected columns from data1
left_join(data1 %>%
select(OBJECTID, year, observed, SMR, hotspot_type),
by = c("OBJECTID", "year"))merged_data$hotspot_type <- factor(merged_data$hotspot_type,
levels = c("High-High", "High-Low", "Low-High", "Low-Low", "Not significant"))
lisa_plot <- ggplot(data = merged_data) +
geom_sf(aes(fill = hotspot_type), color = "black") + # Fill regions based on hotspot type and set borders to black
scale_fill_manual(values = c("High-High" = "red", "Low-Low" = "blue",
"High-Low" = "pink", "Low-High" = "green",
"Not significant" = "grey")) + # Define manual color for each hotspot type
facet_wrap(~ year, ncol = 2) + # Facet by 'year', with two columns
labs(title = "",
fill = "Hotspot Type") +
theme_Publication() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line = element_blank(),
strip.background = element_blank(), # Remove facet label backgrounds
strip.text = element_text(color = "black", size = 12), # Style facet labels
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank(), # Remove y-axis labels
axis.ticks = element_blank(), # Remove axis ticks
legend.title = element_blank(),
legend.position = "right", # Set legend position to top right
legend.justification = "top",
legend.text = element_text(size = 8), # Smaller legend text
legend.key.size = unit(0.5, "cm"),
plot.margin = grid::unit(c(0,0,0,0), "mm")) # Justify the
lisa_plot
merged_data <- merged_data %>%
mutate(
observed_category = cut(observed,
breaks = c(0, 5, 10, Inf),
labels = c("1-5", "6-10", "11+"),
right = TRUE) # intervals are closed on the right
)
# Create a spatial plot with faceting by year, arranged in two columns
oplot <- ggplot(data = merged_data) +
geom_sf(aes(fill = observed_category), color = "black") + # Fill regions based on observed category
scale_fill_manual(values = c("1-5" = "blue", "6-10" = "orange", "11+" = "red"), # Custom colors
name = "Observed Cases", # Rename legend title
guide = guide_legend(title.position = "top")) +
facet_wrap(~ year, ncol = 2) + # Facet by 'year', with two columns
labs(title = "") +
theme_Publication() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line = element_blank(),
strip.background = element_blank(), # Remove facet label backgrounds
strip.text = element_text(color = "black", size = 12), # Style facet labels
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank(), # Remove y-axis labels
axis.ticks = element_blank(), # Remove axis ticks
legend.title = element_blank(),
legend.position = "right", # Set legend position to top right
legend.justification = "top",
legend.text = element_text(size = 8), # Smaller legend text
legend.key.size = unit(0.5, "cm"),
plot.margin = grid::unit(c(0,0,0,0), "mm"))
oplot
Model fitting
Spatial convolution model (1a)
formula1a<-adjusted_observed~1+f(idarea1, model="bym",graph=g)
result1a<-inla(formula1a,family="poisson",data=gdf,
E=E,control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result1a)
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.611, Running = 1.57, Post = 0.298, Total = 2.48
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.003 0.014 -0.03 -0.003 0.025 -0.003 0
Random effects:
Name Model
idarea1 BYM model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 37.85 13.24 18.05 35.81
Precision for idarea1 (spatial component) 17.73 9.76 5.79 15.50
0.975quant mode
Precision for idarea1 (iid component) 69.50 32.05
Precision for idarea1 (spatial component) 42.90 11.91
Deviance Information Criterion (DIC) ...............: 10672.31
Deviance Information Criterion (DIC, saturated) ....: 2639.15
Effective number of parameters .....................: 165.26
Watanabe-Akaike information criterion (WAIC) ...: 10696.48
Effective number of parameters .................: 176.82
Marginal log-Likelihood: -5379.66
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<-adjusted_observed~1+f(idarea1,model="bym",graph=g) + poor + POP_DENS
result1b<-inla(formula1b,family="poisson",data=gdf,
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.363, Running = 1.64, Post = 0.415, Total = 2.42
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.021 0.022 -0.064 -0.021 0.022 -0.021 0
poor 0.005 0.003 -0.001 0.005 0.012 0.005 0
POP_DENS 0.000 0.005 -0.010 0.000 0.009 0.000 0
Random effects:
Name Model
idarea1 BYM model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 37.99 13.92 17.12 35.85
Precision for idarea1 (spatial component) 18.20 10.90 5.63 15.54
0.975quant mode
Precision for idarea1 (iid component) 71.14 31.91
Precision for idarea1 (spatial component) 46.73 11.52
Deviance Information Criterion (DIC) ...............: 10672.61
Deviance Information Criterion (DIC, saturated) ....: 2639.46
Effective number of parameters .....................: 166.19
Watanabe-Akaike information criterion (WAIC) ...: 10697.56
Effective number of parameters .................: 178.29
Marginal log-Likelihood: -5396.29
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 (2a)
formula2a<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+idtime
result2a<-inla(formula2a,family="poisson",data=gdf,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result2a)
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.341, Running = 1.63, Post = 0.404, Total = 2.38
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.004 0.021 -0.037 0.004 0.045 0.004 0
idtime -0.001 0.003 -0.007 -0.001 0.005 -0.001 0
Random effects:
Name Model
idarea1 BYM model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 37.81 13.21 18.05 35.78
Precision for idarea1 (spatial component) 17.75 9.76 5.80 15.52
0.975quant mode
Precision for idarea1 (iid component) 69.36 32.03
Precision for idarea1 (spatial component) 42.92 11.93
Deviance Information Criterion (DIC) ...............: 10674.11
Deviance Information Criterion (DIC, saturated) ....: 2640.95
Effective number of parameters .....................: 166.26
Watanabe-Akaike information criterion (WAIC) ...: 10698.66
Effective number of parameters .................: 178.04
Marginal log-Likelihood: -5388.84
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<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+idtime + poor + POP_DENS
result2b<-inla(formula2b,family="poisson",data=gdf,
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.343, Running = 1.65, Post = 0.423, Total = 2.42
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.087 0.049 -0.183 -0.087 0.010 -0.087 0
idtime 0.008 0.005 -0.002 0.008 0.017 0.008 0
poor 0.012 0.005 0.001 0.012 0.023 0.012 0
POP_DENS 0.000 0.005 -0.009 0.001 0.010 0.001 0
Random effects:
Name Model
idarea1 BYM model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 36.34 12.13 17.65 34.64
Precision for idarea1 (spatial component) 19.73 11.81 6.14 16.84
0.975quant mode
Precision for idarea1 (iid component) 64.82 31.50
Precision for idarea1 (spatial component) 50.68 12.48
Deviance Information Criterion (DIC) ...............: 10673.86
Deviance Information Criterion (DIC, saturated) ....: 2640.70
Effective number of parameters .....................: 167.32
Watanabe-Akaike information criterion (WAIC) ...: 10698.76
Effective number of parameters .................: 179.25
Marginal log-Likelihood: -5403.95
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)
formula3a<-adjusted_observed~1+f(idarea1, model="bym",graph=g) + f(idtime,model="rw2")
result3a<-inla(formula3a,family="poisson",data=gdf,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result3a)
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.362, Running = 1.14, Post = 0.445, Total = 1.95
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.003 0.014 -0.031 -0.003 0.025 -0.003 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 37.58 13.22 17.80 35.54
Precision for idarea1 (spatial component) 18.13 10.39 5.70 15.69
Precision for idtime 30648.41 27432.85 4017.49 22816.65
0.975quant mode
Precision for idarea1 (iid component) 69.15 31.79
Precision for idarea1 (spatial component) 45.08 11.84
Precision for idtime 103272.71 11099.27
Deviance Information Criterion (DIC) ...............: 10675.33
Deviance Information Criterion (DIC, saturated) ....: 2642.17
Effective number of parameters .....................: 167.18
Watanabe-Akaike information criterion (WAIC) ...: 10700.48
Effective number of parameters .................: 179.37
Marginal log-Likelihood: -5386.39
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 2 + cov (3b)
formula3b<-adjusted_observed~1+f(idarea1, model="bym",graph=g) + f(idtime,model="rw2") + poor + POP_DENS
result3b<-inla(formula3b,family="poisson",data=gdf,
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.365, Running = 1.2, Post = 0.51, Total = 2.08
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.050 0.028 -0.105 -0.050 0.006 -0.050 0
poor 0.013 0.006 0.002 0.013 0.024 0.013 0
POP_DENS 0.001 0.005 -0.009 0.001 0.010 0.001 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 36.48 12.84 17.00 34.58
Precision for idarea1 (spatial component) 20.40 13.11 5.82 17.08
Precision for idtime 28667.79 25772.92 3592.58 21279.86
0.975quant mode
Precision for idarea1 (iid component) 66.87 31.08
Precision for idarea1 (spatial component) 54.94 12.23
Precision for idtime 96903.13 10033.18
Deviance Information Criterion (DIC) ...............: 10674.67
Deviance Information Criterion (DIC, saturated) ....: 2641.51
Effective number of parameters .....................: 167.96
Watanabe-Akaike information criterion (WAIC) ...: 10700.19
Effective number of parameters .................: 180.35
Marginal log-Likelihood: -5401.27
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)')Bernadinielli (4a)
formula4a <- adjusted_observed ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid") + idtime
result4a <- inla(formula4a,
family = "poisson", data = gdf, E = E,
control.predictor = list(compute = TRUE),
control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE)
)
summary(result4a)
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.362, Running = 1.08, Post = 0.465, Total = 1.91
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.006 0.021 -0.036 0.006 0.048 0.006 0
idtime -0.002 0.003 -0.008 -0.002 0.004 -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) 38.40 13.72 18.10 36.21
Precision for idarea (spatial component) 18.90 10.90 5.78 16.35
Precision for idarea1 20857.79 33672.03 2587.37 11372.53
0.975quant mode
Precision for idarea (iid component) 71.36 32.21
Precision for idarea (spatial component) 47.10 12.29
Precision for idarea1 100384.40 5139.36
Deviance Information Criterion (DIC) ...............: 10663.45
Deviance Information Criterion (DIC, saturated) ....: 2630.29
Effective number of parameters .....................: 186.24
Watanabe-Akaike information criterion (WAIC) ...: 10687.76
Effective number of parameters .................: 195.33
Marginal log-Likelihood: -5387.71
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 <- adjusted_observed ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid") + idtime + poor + POP_DENS
result4b <- inla(formula4b,
family = "poisson", data = gdf, 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.368, Running = 1.17, Post = 0.607, Total = 2.15
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.084 0.050 -0.182 -0.084 0.014 -0.084 0
idtime 0.007 0.005 -0.003 0.007 0.017 0.007 0
poor 0.012 0.006 0.001 0.012 0.023 0.012 0
POP_DENS 0.000 0.005 -0.009 0.000 0.009 0.000 0
Random effects:
Name Model
idarea BYM model
idarea1 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 37.22 13.19 17.39 35.22
Precision for idarea (spatial component) 21.11 13.46 5.96 17.73
Precision for idarea1 22971.67 39836.67 2625.33 11949.09
0.975quant mode
Precision for idarea (iid component) 68.60 31.54
Precision for idarea (spatial component) 56.46 12.72
Precision for idarea1 114590.26 5203.93
Deviance Information Criterion (DIC) ...............: 10663.60
Deviance Information Criterion (DIC, saturated) ....: 2630.45
Effective number of parameters .....................: 187.06
Watanabe-Akaike information criterion (WAIC) ...: 10688.14
Effective number of parameters .................: 196.18
Marginal log-Likelihood: -5402.82
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
formula5a<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+ f(idtime,model="rw2")+f(id2,model="iid")
result5a<-inla(formula5a,family="poisson",data=gdf,
E=E,control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE))
summary(result5a)
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.362, Running = 2.62, Post = 0.636, Total = 3.62
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.01 0.014 -0.038 -0.01 0.017 -0.01 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 42.64 17.14 18.21 39.63
Precision for idarea1 (spatial component) 17.37 9.97 5.37 15.04
Precision for idtime 29919.55 26567.90 3675.34 22320.29
Precision for id2 60.17 17.50 34.70 57.27
0.975quant mode
Precision for idarea1 (iid component) 84.53 34.20
Precision for idarea1 (spatial component) 43.16 11.33
Precision for idtime 100195.70 10366.18
Precision for id2 102.80 51.33
Deviance Information Criterion (DIC) ...............: 10646.53
Deviance Information Criterion (DIC, saturated) ....: 2613.38
Effective number of parameters .....................: 382.31
Watanabe-Akaike information criterion (WAIC) ...: 10660.50
Effective number of parameters .................: 345.15
Marginal log-Likelihood: -5383.00
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 2 + idd space-time interaction + covariates
formula5b<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+ f(idtime,model="rw2")+f(id2,model="iid") + poor + POP_DENS
result5b<-inla(formula5b,family="poisson",data=gdf,
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.413, Running = 2.78, Post = 0.593, Total = 3.79
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.059 0.029 -0.116 -0.059 -0.002 -0.059 0
poor 0.013 0.006 0.001 0.013 0.025 0.013 0
POP_DENS 0.001 0.005 -0.009 0.001 0.010 0.001 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea1 (iid component) 40.52 15.68 17.71 37.91
Precision for idarea1 (spatial component) 19.64 12.32 5.60 16.59
Precision for idtime 28477.16 25467.02 3700.13 21204.86
Precision for id2 62.05 18.76 35.28 58.80
0.975quant mode
Precision for idarea1 (iid component) 78.44 33.15
Precision for idarea1 (spatial component) 51.93 11.98
Precision for idtime 95893.89 10253.22
Precision for id2 108.17 52.16
Deviance Information Criterion (DIC) ...............: 10644.69
Deviance Information Criterion (DIC, saturated) ....: 2611.53
Effective number of parameters .....................: 380.01
Watanabe-Akaike information criterion (WAIC) ...: 10660.90
Effective number of parameters .................: 345.10
Marginal log-Likelihood: -5398.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)')Analyzing M5b - best fitted
result5b$summary.fixed mean sd 0.025quant 0.5quant 0.975quant
(Intercept) -0.058656493 0.029227773 -0.116164426 -0.058594316 -0.001501804
poor 0.013000441 0.006105956 0.001044326 0.012993348 0.024997089
POP_DENS 0.001099214 0.004811746 -0.008503049 0.001153255 0.010395619
mode kld
(Intercept) -0.058593814 1.208159e-09
poor 0.012993354 3.797174e-10
POP_DENS 0.001153606 3.234118e-08result5b$summary.hyperpar mean sd 0.025quant
Precision for idarea1 (iid component) 40.52161 15.67830 17.708705
Precision for idarea1 (spatial component) 19.64217 12.32440 5.597736
Precision for idtime 28477.15977 25467.02263 3700.131884
Precision for id2 62.05037 18.75987 35.284526
0.5quant 0.975quant mode
Precision for idarea1 (iid component) 37.91351 78.43653 33.14612
Precision for idarea1 (spatial component) 16.58646 51.92822 11.98352
Precision for idtime 21204.85864 95893.88559 10253.21643
Precision for id2 58.80068 108.16767 52.15960summary(result5b$cpo$cpo - result5a$cpo$cpo) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-0.016072 -0.001453 0.000044 -0.000027 0.001605 0.011227 27 plot(result5b$cpo$cpo, result5a$cpo$cpo, xlab = "Model 1 CPO", ylab = "Model 2 CPO", main = "CPO Comparison")
abline(a = 0, b = 1, col = "red")
Sensitivity analysis
formula5b_s2 <- adjusted_observed ~ 1 +
f(idarea1, model = "bym", graph = g, hyper = list(
prec.unstruct = list(prior = "loggamma", param = c(0.01, 0.01)),
prec.spatial = list(prior = "loggamma", param = c(0.01, 0.01))
)) +
f(idtime, model = "rw2", hyper = list(
prec = list(prior = "pc.prec", param = c(0.01, 0.01))
)) +
f(id2, model = "iid", hyper = list(
prec = list(prior = "loggamma", param = c(0.01, 0.01))
)) +
poor + POP_DENS
result5b_s2 <- inla(formula5b_s2, family = "poisson", data = gdf, E = E, control.predictor = list(compute = TRUE), control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE) )
summary(result5b_s2)
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.409, Running = 2.72, Post = 0.714, Total = 3.84
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.056 0.029 -0.112 -0.056 0.000 -0.056 0
poor 0.012 0.006 0.001 0.012 0.024 0.012 0
POP_DENS 0.001 0.005 -0.009 0.001 0.010 0.001 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant
Precision for idarea1 (iid component) 4.04e+01 1.38e+01 19.62
Precision for idarea1 (spatial component) 1.64e+01 8.10e+00 6.03
Precision for idtime 4.90e+06 4.54e+07 18666.84
Precision for id2 5.38e+01 1.36e+01 32.89
0.5quant 0.975quant mode
Precision for idarea1 (iid component) 38.36 7.31e+01 34.57
Precision for idarea1 (spatial component) 14.71 3.71e+01 11.83
Precision for idtime 433031.31 3.19e+07 35817.39
Precision for id2 51.90 8.61e+01 47.88
Deviance Information Criterion (DIC) ...............: 10643.59
Deviance Information Criterion (DIC, saturated) ....: 2610.44
Effective number of parameters .....................: 399.04
Watanabe-Akaike information criterion (WAIC) ...: 10655.67
Effective number of parameters .................: 356.01
Marginal log-Likelihood: -5398.13
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)')formula5b_s3 <- adjusted_observed ~ 1 +
f(idarea1, model = "bym", graph = g, hyper = list(
prec.unstruct = list(prior = "loggamma", param = c(2, 0.5)),
prec.spatial = list(prior = "loggamma", param = c(2, 0.5))
)) +
f(idtime, model = "rw2", hyper = list(
prec = list(prior = "pc.prec", param = c(2, 0.5))
)) +
f(id2, model = "iid", hyper = list(
prec = list(prior = "loggamma", param = c(2, 0.5))
)) +
poor + POP_DENS
result5b_s3 <- inla(formula5b_s3, family = "poisson", data = gdf, E = E, control.predictor = list(compute = TRUE), control.compute=list(dic=TRUE,cpo=TRUE,waic=TRUE, return.marginals.predictor=TRUE) )
summary(result5b_s3)
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.383, Running = 2.7, Post = 1.11, Total = 4.2
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.064 0.032 -0.126 -0.063 -0.001 -0.063 0
poor 0.013 0.007 0.000 0.013 0.025 0.013 0
POP_DENS 0.001 0.005 -0.010 0.001 0.011 0.001 0
Random effects:
Name Model
idarea1 BYM model
idtime RW2 model
id2 IID model
Model hyperparameters:
mean sd 0.025quant
Precision for idarea1 (iid component) 26.57 4.33e+00 19.01
Precision for idarea1 (spatial component) 14.13 3.58e+00 8.33
Precision for idtime 178158.47 1.86e+06 444.48
Precision for id2 32.10 4.22e+00 24.58
0.5quant 0.975quant mode
Precision for idarea1 (iid component) 26.25 3.60e+01 25.63
Precision for idarea1 (spatial component) 13.72 2.23e+01 12.95
Precision for idtime 13481.76 1.14e+06 758.69
Precision for id2 31.84 4.12e+01 31.32
Deviance Information Criterion (DIC) ...............: 10658.78
Deviance Information Criterion (DIC, saturated) ....: 2625.62
Effective number of parameters .....................: 522.13
Watanabe-Akaike information criterion (WAIC) ...: 10638.51
Effective number of parameters .................: 421.35
Marginal log-Likelihood: -5416.42
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)')Marginal posterior distribution
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)
}
# Ensure the 'Model' factor is ordered as desired for the legend
combined_data$Model <- factor(combined_data$Model, levels = model_names)
# Define line types and widths for better differentiation
line_types <- c("solid", "dashed", "dotdash", "longdash")
# Generate the plot
p <- ggplot(combined_data, aes(x, y, color = Model)) +
geom_line(aes(linetype = Model), size = 1.0) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_colour_manual(values = c("#386cb0", "#fdb462", "#7fc97f")) +
scale_linetype_manual(values = line_types) +
theme_Publication(base_size = 12) +
theme(legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 8)) +
labs(title = "",
x = NULL, y = "Probability density")
return(p)
}
# Usage
result_list <- list(result5b, result5b_s2, result5b_s3)
model_names <- c("Default priors, log-gamma (1, 0.0005)", "Best-fitted, log-gamma (0.01, 0.01)", "log-gamma (2, 0.5)")
plot <- 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.print(plot)
More models
Here we explore more model
Space time interaction type II
formula.intII<- adjusted_observed ~ f(idarea,model="bym",graph=g) +
f(idtime,model="rw2") +
f(idtime1,model="iid") +
f(idarea.int,model="iid", group=idyear.int,
control.group=list(model="rw2"))
mod.intII <- inla(formula.intII,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE,cpo=TRUE, waic=TRUE))
summary(mod.intII)
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.446, Running = 59, Post = 0.73, Total = 60.2
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.018 1.429 -2.819 -0.018 2.784 -0.018 0
Random effects:
Name Model
idarea BYM model
idtime RW2 model
idtime1 IID model
idarea.int IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 364.73 8.95 343.77 366.13
Precision for idarea (spatial component) 295.62 7.15 278.87 296.73
Precision for idtime 3741.80 79.88 3555.82 3753.34
Precision for idtime1 30059.44 906.15 28889.12 29881.86
Precision for idarea.int 238.07 6.96 221.75 239.27
0.975quant mode
Precision for idarea (iid component) 377.02 374.03
Precision for idarea (spatial component) 305.46 303.03
Precision for idtime 3854.23 3808.43
Precision for idtime1 32252.77 29028.29
Precision for idarea.int 247.47 245.73
Deviance Information Criterion (DIC) ...............: 10805.22
Deviance Information Criterion (DIC, saturated) ....: 2772.06
Effective number of parameters .....................: 522.00
Watanabe-Akaike information criterion (WAIC) ...: 10807.10
Effective number of parameters .................: 438.32
Marginal log-Likelihood: -6209.97
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)')Space time interaction type III
formula.intIII<- adjusted_observed ~ f(idarea,model="bym",graph=g) +
f(idtime,model="rw2") +
f(idtime1,model="iid") +
f(idyear.int,model="iid", group=idarea.int,
control.group=list(model="besag",
graph=g))
mod.intIII <- inla(formula.intIII,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE,cpo=TRUE, waic=TRUE))
summary(mod.intIII)
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.429, Running = 128, Post = 0.737, Total = 129
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.002 4.676 -9.167 0.002 9.172 0.002 0
Random effects:
Name Model
idarea BYM model
idtime RW2 model
idtime1 IID model
idyear.int IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 49.23 0.402 48.46 49.22
Precision for idarea (spatial component) 10.39 0.169 10.01 10.41
Precision for idtime 314.06 8.222 294.68 315.49
Precision for idtime1 379.98 2.139 376.20 380.06
Precision for idyear.int 65.71 0.636 64.34 65.72
0.975quant mode
Precision for idarea (iid component) 50.06 49.18
Precision for idarea (spatial component) 10.64 10.54
Precision for idtime 324.88 322.22
Precision for idtime1 384.75 379.27
Precision for idyear.int 66.86 65.91
Deviance Information Criterion (DIC) ...............: 10643.65
Deviance Information Criterion (DIC, saturated) ....: 2610.50
Effective number of parameters .....................: 299.30
Watanabe-Akaike information criterion (WAIC) ...: 10665.59
Effective number of parameters .................: 285.93
Marginal log-Likelihood: -6307.06
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)')Space time interaction type III + covariates
formula.intIII<- adjusted_observed ~ f(idarea,model="bym",graph=g) +
f(idtime,model="rw2") +
f(idtime1,model="iid") +
f(idyear.int,model="iid", group=idarea.int,
control.group=list(model="besag",
graph=g)) + poor + POP_DENS
mod.intIII <- inla(formula.intIII,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE,cpo=TRUE, waic=TRUE))
summary(mod.intIII)
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.472, Running = 212, Post = 0.811, Total = 213
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.051 232520.772 -456057.739 -0.051 456057.636 -0.051 0
poor 0.016 0.007 0.001 0.016 0.031 0.016 0
POP_DENS 0.002 0.004 -0.006 0.002 0.011 0.002 0
Random effects:
Name Model
idarea BYM model
idtime RW2 model
idtime1 IID model
idyear.int IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for idarea (iid component) 29.31 0.254 28.76 29.31
Precision for idarea (spatial component) 29.91 0.260 29.43 29.91
Precision for idtime 133.85 1.708 129.96 134.00
Precision for idtime1 75.40 0.738 74.12 75.37
Precision for idyear.int 88.41 1.595 86.17 88.19
0.975quant mode
Precision for idarea (iid component) 29.76 29.40
Precision for idarea (spatial component) 30.47 29.85
Precision for idtime 136.52 135.02
Precision for idtime1 77.03 75.10
Precision for idyear.int 92.19 87.01
Deviance Information Criterion (DIC) ...............: 10649.66
Deviance Information Criterion (DIC, saturated) ....: 2616.50
Effective number of parameters .....................: 274.19
Watanabe-Akaike information criterion (WAIC) ...: 10674.16
Effective number of parameters .................: 268.06
Marginal log-Likelihood: -6324.84
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)')Space time interaction type IV
formula.intIV<- adjusted_observed ~ 1+ f(idarea,model="bym",graph=g) +
f(idtime,model="rw2") +
f(idtime1,model="iid") +
f(idarea.int,model="besag", graph=g,
group=idyear.int,
control.group=list(model="rw2"))
mod.intIV <- inla(formula.intIV,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE,cpo=TRUE, waic=TRUE))