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 84Calculate 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’s I Calculation
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.452, Running = 1.3, Post = 0.0919, Total = 1.84
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.307, Running = 1.3, Post = 0.0878, Total = 1.7
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 + temporal random walk 1 (2a)
formula2a<-adjusted_observed~1+f(idarea1, model="bym",graph=g) + f(idtime,model="rw2")
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.308, Running = 0.851, Post = 0.0866, Total = 1.25
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.14 10.39 5.70 15.69
Precision for idtime 30648.58 27433.81 4017.31 22816.45
0.975quant mode
Precision for idarea1 (iid component) 69.15 31.79
Precision for idarea1 (spatial component) 45.08 11.84
Precision for idtime 103275.52 11098.80
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 + covariates (2b)
formula2b<-adjusted_observed~1+f(idarea1, model="bym",graph=g) + f(idtime,model="rw2") + 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.32, Running = 0.868, Post = 0.0947, Total = 1.28
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.74 25772.68 3592.63 21279.92
0.975quant mode
Precision for idarea1 (iid component) 66.87 31.08
Precision for idarea1 (spatial component) 54.94 12.23
Precision for idtime 96902.41 10033.33
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)')Space-time interaction additive model
formula3a<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+ f(idtime,model="rw2")+f(id2,model="iid")
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.33, Running = 1.74, Post = 0.191, Total = 2.26
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.15 18.21 39.63
Precision for idarea1 (spatial component) 17.37 9.97 5.37 15.04
Precision for idtime 29919.68 26569.62 3674.70 22319.68
Precision for id2 60.17 17.50 34.70 57.27
0.975quant mode
Precision for idarea1 (iid component) 84.53 34.21
Precision for idarea1 (spatial component) 43.16 11.33
Precision for idtime 100200.56 10364.56
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)')Space-time interaction additive model + covariates
formula3b<-adjusted_observed~1+f(idarea1,model="bym",graph=g)+ f(idtime,model="rw2")+f(id2,model="iid") + 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.327, Running = 1.89, Post = 0.172, Total = 2.39
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.53 15.69 17.70 37.91
Precision for idarea1 (spatial component) 19.64 12.34 5.59 16.59
Precision for idtime 28466.17 25427.95 3705.47 21208.31
Precision for id2 62.05 18.76 35.29 58.80
0.975quant mode
Precision for idarea1 (iid component) 78.48 33.14
Precision for idarea1 (spatial component) 51.96 11.98
Precision for idtime 95776.48 10267.55
Precision for id2 108.16 52.16
Deviance Information Criterion (DIC) ...............: 10644.68
Deviance Information Criterion (DIC, saturated) ....: 2611.52
Effective number of parameters .....................: 380.02
Watanabe-Akaike information criterion (WAIC) ...: 10660.90
Effective number of parameters .................: 345.11
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)')Space time interaction type II
formula4a<- 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"))
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.354, Running = 99.3, Post = 0.222, Total = 99.9
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.01 97080.73 -190402.5 -0.01 190402.5 -0.01 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) 1982.05 10.97 1964.00 1982.72
Precision for idarea (spatial component) 2223.54 13.74 2202.24 2224.30
Precision for idtime 20624.21 79.43 20449.88 20623.47
Precision for idtime1 24164.70 347.46 23406.60 24181.37
Precision for idarea.int 296.20 4.42 289.78 295.68
0.975quant mode
Precision for idarea (iid component) 2007.17 1977.19
Precision for idarea (spatial component) 2256.01 2215.36
Precision for idtime 20785.33 20628.90
Precision for idtime1 24766.36 24312.26
Precision for idarea.int 306.60 292.66
Deviance Information Criterion (DIC) ...............: 10803.28
Deviance Information Criterion (DIC, saturated) ....: 2770.12
Effective number of parameters .....................: 511.07
Watanabe-Akaike information criterion (WAIC) ...: 10808.15
Effective number of parameters .................: 432.63
Marginal log-Likelihood: -6206.28
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 II + covariates
formula4b<- 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")) + + 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))
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.361, Running = 41, Post = 0.178, Total = 41.5
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.039 142215.063 -283704.424 0.040 283704.505 0.041 0
poor 0.009 0.013 -0.016 0.009 0.035 0.009 0
POP_DENS -0.050 0.031 -0.110 -0.050 0.010 -0.050 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) 48.16 0.589 47.34 48.10
Precision for idarea (spatial component) 59.07 0.408 58.27 59.07
Precision for idtime 58.70 0.753 57.27 58.68
Precision for idtime1 61.87 1.740 57.91 62.07
Precision for idarea.int 86.50 1.268 83.78 86.58
0.975quant mode
Precision for idarea (iid component) 49.57 47.64
Precision for idarea (spatial component) 59.92 59.03
Precision for idtime 60.27 58.59
Precision for idtime1 64.54 63.23
Precision for idarea.int 88.72 86.99
Deviance Information Criterion (DIC) ...............: 10837.07
Deviance Information Criterion (DIC, saturated) ....: 2803.91
Effective number of parameters .....................: 601.75
Watanabe-Akaike information criterion (WAIC) ...: 10820.21
Effective number of parameters .................: 481.54
Marginal log-Likelihood: -6254.34
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
formula5a<- 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))
result5a <- inla(formula5a,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
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.44, Running = 47.2, Post = 0.168, Total = 47.8
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.004 207977.1 -419355.9 -0.01 419355.9 -0.01 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) 34.32 0.718 32.65 34.41
Precision for idarea (spatial component) 14.73 0.109 14.52 14.73
Precision for idtime 168.67 1.604 165.08 168.70
Precision for idtime1 207.49 2.781 201.02 207.77
Precision for idyear.int 75.62 0.588 74.39 75.63
0.975quant mode
Precision for idarea (iid component) 35.34 34.91
Precision for idarea (spatial component) 14.96 14.71
Precision for idtime 171.37 169.45
Precision for idtime1 211.57 209.74
Precision for idyear.int 76.74 75.73
Deviance Information Criterion (DIC) ...............: 10644.39
Deviance Information Criterion (DIC, saturated) ....: 2611.23
Effective number of parameters .....................: 290.04
Watanabe-Akaike information criterion (WAIC) ...: 10667.24
Effective number of parameters .................: 279.37
Marginal log-Likelihood: -6307.86
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
formula5b<- 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
result5b <- inla(formula5b,family="poisson",data=gdf,E=E,
control.predictor=list(compute=TRUE),
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.402, Running = 152, Post = 0.195, Total = 152
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.048 183366.403 -359582.005 -0.048 359581.908 -0.048 0
poor 0.016 0.008 0.001 0.016 0.030 0.016 0
POP_DENS 0.002 0.005 -0.007 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) 31.32 0.813 29.51 31.37
Precision for idarea (spatial component) 20.98 0.984 19.52 20.85
Precision for idtime 515.28 11.566 488.32 516.85
Precision for idtime1 614.27 6.412 602.76 614.27
Precision for idyear.int 85.15 1.106 83.53 85.11
0.975quant mode
Precision for idarea (iid component) 32.66 31.78
Precision for idarea (spatial component) 23.30 20.25
Precision for idtime 531.57 524.95
Precision for idtime1 628.44 612.15
Precision for idyear.int 87.78 84.32
Deviance Information Criterion (DIC) ...............: 10649.11
Deviance Information Criterion (DIC, saturated) ....: 2615.95
Effective number of parameters .....................: 273.49
Watanabe-Akaike information criterion (WAIC) ...: 10673.56
Effective number of parameters .................: 267.47
Marginal log-Likelihood: -6319.62
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)')Comparing model performance
# Apply inla.group.cv to each result
apply_inla_group_cv <- function(result) {
result$loocv <- inla.group.cv(result, num.level.sets = -1)
result$logcv.m3 <- inla.group.cv(result, num.level.sets = 3)
result$logcv.m5 <- inla.group.cv(result, num.level.sets = 5)
result$logcv.m10 <- inla.group.cv(result, num.level.sets = 10)
return(result)
}
# Apply the function to each result in the list
results <- list(result1a, result1b, result2a, result2b, result3a, result3b, result4a, result4b, result5a, result5b)
results <- lapply(results, apply_inla_group_cv)
names(results) <- c("Model 1a", "Model 1b", "Model 2a", "Model 2b", "Model 3a", "Model 3b", "Model 4a", "Model 4b", "Model 5a", "Model 5b")
# Define functions for extracting metrics
DIC <- function(x) {
data.frame(mean.deviance = x$dic$mean.deviance,
p.eff.DIC = x$dic$p.eff,
DIC = x$dic$dic,
WAIC = x$waic$waic,
p.WAIC = x$waic$p.eff)
}
compute_loocv <- function(x) {
cv_values <- x$loocv$cv
cv_values <- cv_values[!is.na(cv_values)] # Remove NA values
if (length(cv_values) > 0) {
-mean(log(cv_values))
} else {
NA
}
}
logcv <- function(cv_values) {
cv_values <- cv_values[!is.na(cv_values)] # Remove NA values
if (length(cv_values) > 0) {
-mean(log(cv_values))
} else {
NA # Handle cases where all values are NA
}
}
compute_logcv <- function(x) {
# Extract logcv values for each num.level.sets
logcv_m3 <- logcv(x$logcv.m3$cv)
logcv_m5 <- logcv(x$logcv.m5$cv)
logcv_m10 <- logcv(x$logcv.m10$cv)
#rounding
logcv_m3_rounded <- if (!is.na(logcv_m3)) round(logcv_m3, 3) else NA
logcv_m5_rounded <- if (!is.na(logcv_m5)) round(logcv_m5, 3) else NA
logcv_m10_rounded <- if (!is.na(logcv_m10)) round(logcv_m10, 3) else NA
list(
logcv_m3 = logcv_m3_rounded,
logcv_m5 = logcv_m5_rounded,
logcv_m10 = logcv_m10_rounded
)
}
extract_group_indices <- function(x, num.level.sets) {
lapply(x[[paste0("logcv.m", num.level.sets)]]$groups, function(group) group$idx)
}
tables <- list()
for (name in names(results)) {
result <- results[[name]]
# Extract DIC and WAIC
dic_waic <- DIC(result)
# print(paste("DIC for", name, ":", dic_waic))
# Compute LOOCV
loocv <- round(compute_loocv(result), 3)
# print(paste("LOOCV for", name, ":", loocv))
# Compute logcv
logcv_results <- compute_logcv(result)
# print(paste("logcv_m3 for", name, ":", logcv_results$logcv_m3))
# print(paste("logcv_m5 for", name, ":", logcv_results$logcv_m5))
# print(paste("logcv_m10 for", name, ":", logcv_results$logcv_m10))
# Extract group indices
# groups_m3 <- extract_group_indices(result, 3)
# groups_m5 <- extract_group_indices(result, 5)
# groups_m10 <- extract_group_indices(result, 10)
# print(paste("Groups_m3 for", name, ":", paste(unlist(groups_m3), collapse = ", ")))
# print(paste("Groups_m5 for", name, ":", paste(unlist(groups_m5), collapse = ", ")))
# print(paste("Groups_m10 for", name, ":", paste(unlist(groups_m10), collapse = ", ")))
# Combine all data into a single table
table <- data.frame(
mean.deviance = dic_waic$mean.deviance,
DIC = dic_waic$DIC,
p.DIC = dic_waic$p.eff,
WAIC = dic_waic$WAIC,
p.WAIC = dic_waic$p.WAIC,
LOOCV = loocv,
logcv_m3 = I(list(logcv_results$logcv_m3)),
logcv_m5 = I(list(logcv_results$logcv_m5)),
logcv_m10 = I(list(logcv_results$logcv_m10))
)
tables[[name]] <- table
}
# Combine all tables into a single data frame
final_table <- do.call(rbind, tables)
# Print the final table
kable(final_table, format = "html", caption = "Model Evaluation Results")| mean.deviance | DIC | p.DIC | WAIC | p.WAIC | LOOCV | logcv_m3 | logcv_m5 | logcv_m10 | |
|---|---|---|---|---|---|---|---|---|---|
| Model 1a | 10507.05 | 10672.31 | 165.2574 | 10696.48 | 176.8209 | 2.395 | 2.463 | 2.463 | 2.463 |
| Model 1b | 10506.42 | 10672.61 | 166.1913 | 10697.56 | 178.2874 | 2.395 | 2.414 | 2.432 | 2.467 |
| Model 2a | 10508.15 | 10675.33 | 167.1821 | 10700.48 | 179.3718 | 2.396 | 2.418 | 2.436 | 2.467 |
| Model 2b | 10506.71 | 10674.67 | 167.9594 | 10700.19 | 180.3550 | 2.396 | 2.419 | 2.443 | 2.467 |
| Model 3a | 10264.22 | 10646.53 | 382.3128 | 10660.50 | 345.1542 | 2.391 | 2.437 | 2.454 | 2.449 |
| Model 3b | 10264.66 | 10644.68 | 380.0167 | 10660.90 | 345.1061 | 2.391 | 2.437 | 2.456 | 2.448 |
| Model 4a | 10292.21 | 10803.28 | 511.0710 | 10808.15 | 432.6329 | 2.431 | 2.561 | 2.844 | 4.299 |
| Model 4b | 10235.32 | 10837.07 | 601.7524 | 10820.21 | 481.5438 | 2.440 | 2.576 | 2.82 | 3.646 |
| Model 5a | 10354.35 | 10644.38 | 290.0367 | 10667.24 | 279.3690 | 2.391 | 2.428 | 2.454 | 2.465 |
| Model 5b | 10375.62 | 10649.11 | 273.4914 | 10673.56 | 267.4711 | 2.392 | 2.431 | 2.458 | 2.466 |
Analyzing M3b - best fitted
result3b$summary.fixed mean sd 0.025quant 0.5quant 0.975quant
(Intercept) -0.058655807 0.029227316 -0.116162673 -0.058593674 -0.001501919
poor 0.013000225 0.006105871 0.001044241 0.012993145 0.024996666
POP_DENS 0.001099171 0.004811779 -0.008503182 0.001153227 0.010395601
mode kld
(Intercept) -0.058593169 1.205954e-09
poor 0.012993152 3.784022e-10
POP_DENS 0.001153583 3.234859e-08result3b$summary.hyperpar mean sd 0.025quant
Precision for idarea1 (iid component) 40.52628 15.69245 17.700308
Precision for idarea1 (spatial component) 19.64479 12.33589 5.592538
Precision for idtime 28466.16957 25427.94895 3705.471339
Precision for id2 62.04805 18.75721 35.287208
0.5quant 0.975quant mode
Precision for idarea1 (iid component) 37.91362 78.48133 33.13852
Precision for idarea1 (spatial component) 16.58510 51.96276 11.97744
Precision for idtime 21208.31244 95776.47907 10267.55074
Precision for id2 58.79856 108.15967 52.15702summary(result3b$cpo$cpo - result3a$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(result3b$cpo$cpo, result3a$cpo$cpo, xlab = "Model 3 CPO", ylab = "Model 2 CPO", main = "CPO Comparison")
abline(a = 0, b = 1, col = "red")
Sensitivity analysis
formula3b_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
result3b_s2 <- inla(formula3b_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(result3b_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.33, Running = 1.83, Post = 0.23, Total = 2.39
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 18667.78
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 433005.09 3.19e+07 35820.57
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)')formula3b_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
result3b_s3 <- inla(formula3b_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(result3b_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.375, Running = 1.8, Post = 0.316, Total = 2.49
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.57e+00 8.33
Precision for idtime 177535.48 1.85e+06 445.10
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 13475.67 1.13e+06 760.35
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.14
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(result3b, result3b_s2, result3b_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)
Space time interaction type IV