Lower Respiratory Tract Infections in England: A Bayesian spatio-temporal analysis
2023-07-05
library(magrittr)
library(utils)
library(ggplot2)
library(sf)
library(sp)
library(rgdal)
library(dplyr)
library(viridis)
library(RColorBrewer)
library(INLA)load("data.RData")
load("map.RData")Descriptive statistics
table1 <- t(sapply(split.default(data[, c(6, 7, 9, 10, 11, 12, 13, 14)], names(data[, c(6, 7, 9, 10, 11, 12, 13, 14)])), function(x) {
x1 <- unlist(x)
data.frame(Mean = mean(x1, na.rm = TRUE), "Standard deviations" = sd(x1, na.rm = TRUE))
})) %>% data.frame()
desired_order <- c(6, 8, 3, 1, 4, 5, 7, 2)
table1 <- table1[desired_order, ]
rownames(table1)[1] <- "Number of hospitalizations"
rownames(table1)[2] <- "Population"
rownames(table1)[3] <- "log(Population)"
rownames(table1)[4] <- "Total greenspace (ha)"
rownames(table1)[5] <- "NDVI"
rownames(table1)[6] <- "NO$_{2}$"
rownames(table1)[7] <- "PM$_{2.5}$"
# convert t in latex
library(xtable)
print(xtable(table1), type = "latex")## % latex table generated in R 4.3.1 by xtable 1.8-4 package
## % Fri Aug 25 04:04:49 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrr}
## \hline
## & Mean & Standard.deviations \\
## \hline
## Number of hospitalizations & 159.86 & 137.30 \\
## Population & 39294.26 & 28557.33 \\
## log(Population) & 10.39 & 0.59 \\
## Total greenspace (ha) & 949.71 & 663.84 \\
## NDVI & 0.54 & 0.07 \\
## NO\$\_\{2\}\$ & 13.64 & 6.58 \\
## PM\$\_\{2.5\}\$ & 9.65 & 1.84 \\
## IMD & 19.64 & 7.98 \\
## \hline
## \end{tabular}
## \end{table}Incidence map
Spatial distribution of child hospital admission rates from 2013 to 2019.
# merge data frames
incidence_map <- merge(map, data[c("Level", "Year", "O")], by.x = "lad19cd", by.y = "Level")
# create bins for fill colour
incidence_map$O.cut <- cut(
incidence_map$O,
breaks = c(0, 75, 200, 350, 500, max(incidence_map$O[!is.na(incidence_map$O)])),
)
pO <- ggplot() +
geom_sf(data = incidence_map, aes(fill = O.cut)) +
theme_bw() +
scale_fill_brewer(name = "Rate of child hospitalization", palette = "Reds", labels = c("<75", "(75, 200]", "(200, 350]", "(350, 500]", "(500, 1465]")) +
guides(fill = guide_legend(title = NULL)) +
theme(
text = element_text(size = 20),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.9,0.2)
) +
facet_wrap(~Year, ncol = 4) +
labs("")
print(pO)
ggsave("Obs_map.png", pO)
Covariates
png("hist_pm25.png")
hist(data$PM25, main=expression(paste("Histogram of ", PM[2.5])),xlab=expression(PM[2.5] (µg/m^3)))
dev.off()## quartz_off_screen
## 2png("hist_no2.png")
hist(data$NO2, main=expression(paste("Histogram of ", NO^2)), xlab=expression(NO^2* (µg/m^3)))
dev.off()## quartz_off_screen
## 2png("hist_greenspace.png")
hist(data$Greenspace, main = expression("Histogram of greenspace"), xlab="Greenspace (ha)")
dev.off()## quartz_off_screen
## 2png("hist_ndvi.png")
hist(data$NDVI, main = expression("Histogram of NDVI"))
dev.off()## quartz_off_screen
## 2png("hist_imd.png")
hist(data$IMD, main = expression("Histogram of IMD"), xlab="IMD")
dev.off()## quartz_off_screen
## 2png("hist_pop.png")
hist(data$Population, main = expression("Histogram of Population"))
dev.off## function (which = dev.cur())
## {
## if (which == 1)
## stop("cannot shut down device 1 (the null device)")
## .External(C_devoff, as.integer(which))
## dev.cur()
## }
## <bytecode: 0x14632c020>
## <environment: namespace:grDevices># Assuming 'data' is your DataFrame
variables <- c("O", "Population", "PM25", "NO2", "Greenspace", "NDVI", "IMD")
# Create a matrix to store mean and standard deviation
summary_matrix <- matrix(0, nrow = length(variables), ncol = 1)
colnames(summary_matrix) <- c("Mean")
# Fill the matrix with mean values
for (i in 1:length(variables)) {
variable <- variables[i]
summary_matrix[i, "Mean"] <- mean(data[[variable]], na.rm = TRUE)
}
# Convert the matrix to a data frame for better presentation
summary_table <- data.frame(Variable = variables, Mean = summary_matrix)
# Calculate and add standard deviation values to the summary table
summary_table <- summary_table %>%
mutate(StandardDeviation = sapply(variables, function(var) {
sd(data[[var]], na.rm = TRUE)
}))
# Format the numbers to avoid scientific notation
summary_table <- summary_table %>%
mutate(across(c("Mean", "StandardDeviation"), function(x) format(x, scientific = FALSE)))
# Print the summary table
print(summary_table)## Variable Mean StandardDeviation
## 1 O 159.8557866 137.30199097
## 2 Population 39294.2617541 28557.33326631
## 3 PM25 9.6455429 1.83884670
## 4 NO2 13.6420834 6.57719763
## 5 Greenspace 949.7100843 663.84343065
## 6 NDVI 0.5403185 0.06982153
## 7 IMD 19.6408133 7.97525676corr<-cor(data$Greenspace, data$NDVI, use = "complete.obs") %>% round(.,2) %>% as.character()
# correlation
corr_gs<-ggplot(data, aes(x = data$Greenspace, y = data$NDVI)) +
geom_point() +
labs(x = "Greenspace", y = "NDVI", title = "Correlation of greenspace measures") +
theme_minimal() +
geom_text(aes(x = 2000, y = 0.35, label = paste("Correlation: ",corr)), color = "blue")
ggsave("corr_gs.png", corr_gs)
corr2<-cor(data$PM25, data$NO2, use = "complete.obs") %>% round(.,2) %>% as.character()
# correlation
corr_env<-ggplot(data, aes(x = data$PM25, y = data$NO2)) +
geom_point() +
labs(x = expression(PM[2.5]), y = expression(NO^2), title = "Correlation of air pollution measures") +
theme_minimal() +
geom_text(aes(x = 7, y = 40, label = paste("Correlation: ",corr2)), color = "blue")
ggsave("corr_env.png", corr_env)
print(corr_gs)
print(corr_env)Morans Test
Model A: Fixed effects model
Determine best combination:
# aggregate data
cases_agg <- data %>%
group_by(Level) %>%
summarise(Observed = sum(O), Expected = sum(E), Population = sum(Population), meanPM25_scaled = mean(PM25_scaled), meanNO2_scaled = mean(NO2_scaled), meanGreenspace_log_scaled = mean(Greenspace_log_scaled), meanNDVI_scaled = mean(NDVI_scaled), meanIMD_scaled = mean(IMD_scaled)) %>%
mutate(log.Population = log(Population)) %>%
dplyr::rename(O = Observed, E = Expected)formula_fixed <- O ~ meanPM25_scaled + meanNDVI_scaled + meanIMD_scaled
fixed.model <- inla(formula_fixed, family = "poisson", data = cases_agg, offset = log.Population, control.compute = list(waic = TRUE))
summary(fixed.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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.679, Running = 0.171, Post = 0.00974, Total = 0.86
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.526 0.002 -5.530 -5.526 -5.523 -5.526 0
## meanPM25_scaled -0.156 0.003 -0.162 -0.156 -0.151 -0.156 0
## meanNDVI_scaled 0.036 0.003 0.030 0.036 0.042 0.036 0
## meanIMD_scaled 0.101 0.002 0.097 0.101 0.105 0.101 0
##
## Watanabe-Akaike information criterion (WAIC) ...: 110941.04
## Effective number of parameters .................: 49332.93
##
## Marginal log-Likelihood: -13823.50
## 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)')formula_fixed <- O ~ meanNO2_scaled + meanNDVI_scaled + meanIMD_scaled
fixed.model <- inla(formula_fixed, family = "poisson", data = cases_agg, offset = log.Population, control.compute = list(waic = TRUE))
summary(fixed.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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.569, Running = 0.171, Post = 0.00498, Total = 0.744
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.521 0.002 -5.524 -5.521 -5.517 -5.521 0
## meanNO2_scaled -0.141 0.003 -0.147 -0.141 -0.136 -0.141 0
## meanNDVI_scaled 0.040 0.003 0.034 0.040 0.047 0.040 0
## meanIMD_scaled 0.150 0.002 0.146 0.150 0.154 0.150 0
##
## Watanabe-Akaike information criterion (WAIC) ...: 212409.74
## Effective number of parameters .................: 99872.94
##
## Marginal log-Likelihood: -14187.68
## 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)')formula_fixed <- O ~ meanPM25_scaled + meanNDVI_scaled + meanIMD_scaled
fixed.model <- inla(formula_fixed, family = "poisson", data = cases_agg, offset = log.Population, control.compute = list(waic = TRUE))
summary(fixed.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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.548, Running = 0.17, Post = 0.00504, Total = 0.723
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.526 0.002 -5.530 -5.526 -5.523 -5.526 0
## meanPM25_scaled -0.156 0.003 -0.162 -0.156 -0.151 -0.156 0
## meanNDVI_scaled 0.036 0.003 0.030 0.036 0.042 0.036 0
## meanIMD_scaled 0.101 0.002 0.097 0.101 0.105 0.101 0
##
## Watanabe-Akaike information criterion (WAIC) ...: 110941.04
## Effective number of parameters .................: 49332.93
##
## Marginal log-Likelihood: -13823.50
## 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)')formula_fixed <- O ~ meanPM25_scaled + meanGreenspace_log_scaled + meanIMD_scaled
fixed.model <- inla(formula_fixed, family = "poisson", data = cases_agg, offset = log.Population, control.compute = list(waic = TRUE))
summary(fixed.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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.548, Running = 0.168, Post = 0.00521, Total = 0.721
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -5.515 0.002 -5.518 -5.515 -5.511 -5.515
## meanPM25_scaled -0.193 0.002 -0.197 -0.193 -0.189 -0.193
## meanGreenspace_log_scaled -0.038 0.002 -0.042 -0.038 -0.035 -0.038
## meanIMD_scaled 0.083 0.002 0.080 0.083 0.086 0.083
## kld
## (Intercept) 0
## meanPM25_scaled 0
## meanGreenspace_log_scaled 0
## meanIMD_scaled 0
##
## Watanabe-Akaike information criterion (WAIC) ...: 95858.05
## Effective number of parameters .................: 41864.03
##
## Marginal log-Likelihood: -13661.00
## 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)')Model B: Spatial + Fixed effects model
ID <- seq(1, 316)
formula_BYM2 <- O ~ meanPM25_scaled + meanGreenspace_log_scaled + meanIMD_scaled + f(ID, model = "bym2", graph = map.adj, hyper = list(
prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)
)
))
sBYM.model <- inla(formula = formula_BYM2, family = "poisson", data = cases_agg, offset = log.Population, control.compute = list(waic = TRUE))
# Obtain the posterior summary statistics
summary(sBYM.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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 = 4.5, Running = 0.249, Post = 0.016, Total = 4.77
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -5.559 0.005 -5.570 -5.559 -5.548 -5.559
## meanPM25_scaled -0.055 0.049 -0.150 -0.055 0.041 -0.055
## meanGreenspace_log_scaled -0.026 0.014 -0.053 -0.026 0.002 -0.026
## meanIMD_scaled 0.106 0.017 0.073 0.106 0.139 0.106
## kld
## (Intercept) 0
## meanPM25_scaled 0
## meanGreenspace_log_scaled 0
## meanIMD_scaled 0
##
## Random effects:
## Name Model
## ID BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 10.439 1.179 8.186 10.415 12.812 9.131
## Phi for ID 0.919 0.051 0.793 0.931 0.986 0.914
##
## Watanabe-Akaike information criterion (WAIC) ...: 3260.03
## Effective number of parameters .................: 158.58
##
## Marginal log-Likelihood: -1986.99
## 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)')Model C: Spatio-temporal and fixed effects
DATA_ST_FE <- merge(map, data[c("Level", "Year", "O", "E", "log.Population", "PM25_scaled", "Greenspace_log_scaled", "IMD_scaled")], by.x = "lad19cd", by.y = "Level")
# Create the ID for year (time)
DATA_ST_FE$ID.time <- DATA_ST_FE$Year - 2012
# Create the ID for space
DATA_ST_FE$ID.space <- rep(seq(1, 316), each = 7)
formula_ST_FE <- O ~ PM25_scaled + Greenspace_log_scaled + IMD_scaled + f(ID.space,
model = "bym2", graph = map.adj,
hyper = list(
prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)
)
)
) +
f(ID.time, model = "rw1", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
)))
st_FE.BYM.model <- inla(
formula = formula_ST_FE, family = "poisson", data = DATA_ST_FE, offset = log.Population,
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(st_FE.BYM.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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 = 4.69, Running = 0.364, Post = 0.0192, Total = 5.08
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.564 0.005 -5.574 -5.564 -5.554 -5.564 0
## PM25_scaled 0.050 0.007 0.037 0.050 0.063 0.050 0
## Greenspace_log_scaled -0.018 0.014 -0.045 -0.018 0.009 -0.018 0
## IMD_scaled 0.089 0.015 0.060 0.089 0.119 0.089 0
##
## Random effects:
## Name Model
## ID.space BYM2 model
## ID.time RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID.space 9.940 1.020 7.984 9.921 11.993 8.950
## Phi for ID.space 0.943 0.041 0.841 0.953 0.992 0.942
## Precision for ID.time 149.466 94.269 34.435 127.470 390.463 6.414
##
## Deviance Information Criterion (DIC) ...............: 22366.66
## Deviance Information Criterion (DIC, saturated) ....: 7738.12
## Effective number of parameters .....................: -1405.12
##
## Watanabe-Akaike information criterion (WAIC) ...: 41839.13
## Effective number of parameters .................: 9243.57
##
## Marginal log-Likelihood: -13234.29
## 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)')Model D: Spatio-temporal and fixed effects with space:time interaction
DATA_ST_FE <- merge(map, data[c("Level", "Year", "O", "E", "log.Population", "PM25_scaled", "Greenspace_log_scaled", "IMD_scaled")], by.x = "lad19cd", by.y = "Level")
# Create the ID for year (time)
DATA_ST_FE$ID.time <- DATA_ST_FE$Year - 2012
# Create the ID for space
DATA_ST_FE$ID.space <- rep(seq(1, 316), each = 7)
DATA_ST_FE$ID.space.time <- seq(1, dim(DATA_ST_FE)[1])
formula_ST_FE_intI <- O ~ PM25_scaled + Greenspace_log_scaled + IMD_scaled + f(ID.space,
model = "bym2", graph = map.adj,
hyper = list(
prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)
)
)
) +
f(ID.time, model = "rw1", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
))) +
f(ID.space.time, model = "iid", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
)))
st_FE_IntI.BYM.model <- inla(
formula = formula_ST_FE_intI, family = "poisson", data = DATA_ST_FE, offset = log.Population,
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(st_FE_IntI.BYM.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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 = 4.73, Running = 0.652, Post = 0.0348, Total = 5.41
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.578 0.006 -5.589 -5.578 -5.567 -5.578 0
## PM25_scaled 0.033 0.016 0.001 0.033 0.065 0.033 0
## Greenspace_log_scaled -0.018 0.014 -0.046 -0.018 0.009 -0.018 0
## IMD_scaled 0.092 0.016 0.061 0.092 0.122 0.092 0
##
## Random effects:
## Name Model
## ID.space BYM2 model
## ID.time RW1 model
## ID.space.time IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for ID.space 10.314 1.046 8.366 10.276 12.480
## Phi for ID.space 0.959 0.041 0.845 0.972 0.996
## Precision for ID.time 167.834 118.080 36.889 137.638 477.136
## Precision for ID.space.time 31.676 1.352 29.052 31.663 34.375
## mode
## Precision for ID.space 9.307
## Phi for ID.space 0.968
## Precision for ID.time 6.581
## Precision for ID.space.time 29.913
##
## Deviance Information Criterion (DIC) ...............: 18784.61
## Deviance Information Criterion (DIC, saturated) ....: 4077.64
## Effective number of parameters .....................: 1780.57
##
## Watanabe-Akaike information criterion (WAIC) ...: 18512.07
## Effective number of parameters .................: 1098.39
##
## Marginal log-Likelihood: -10537.35
## 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)')Model E: Mixed effects model w/ interaction + fixed interations
formula_ST_intI_FE_int <- O ~ PM25_scaled + Greenspace_log_scaled + IMD_scaled + PM25_scaled:Greenspace_log_scaled +
f(ID.space,
model = "bym2", graph = map.adj,
hyper = list(
prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)
)
)
) +
f(ID.time, model = "rw1", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
))) +
f(ID.space.time, model = "iid", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
)))
stIntI.FE.int.BYM.model <- inla(
formula = formula_ST_intI_FE_int, family = "poisson", data = DATA_ST_FE, offset = log.Population,
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(stIntI.FE.int.BYM.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))##
## 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 = 4.95, Running = 0.894, Post = 0.0322, Total = 5.88
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -5.578 0.005 -5.589 -5.578 -5.567
## PM25_scaled 0.031 0.016 -0.001 0.031 0.063
## Greenspace_log_scaled -0.019 0.014 -0.046 -0.019 0.009
## IMD_scaled 0.092 0.016 0.062 0.092 0.123
## PM25_scaled:Greenspace_log_scaled -0.001 0.007 -0.015 -0.001 0.012
## mode kld
## (Intercept) -5.578 0
## PM25_scaled 0.031 0
## Greenspace_log_scaled -0.019 0
## IMD_scaled 0.092 0
## PM25_scaled:Greenspace_log_scaled -0.001 0
##
## Random effects:
## Name Model
## ID.space BYM2 model
## ID.time RW1 model
## ID.space.time IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for ID.space 9.975 0.966 8.105 9.963 11.90
## Phi for ID.space 0.976 0.048 0.843 0.992 1.00
## Precision for ID.time 140.420 111.609 10.730 110.670 414.54
## Precision for ID.space.time 31.665 1.362 29.126 31.616 34.49
## mode
## Precision for ID.space 9.075
## Phi for ID.space 0.991
## Precision for ID.time 2.335
## Precision for ID.space.time 29.998
##
## Deviance Information Criterion (DIC) ...............: 18784.31
## Deviance Information Criterion (DIC, saturated) ....: 4077.34
## Effective number of parameters .....................: 1781.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 18509.92
## Effective number of parameters .................: 1097.69
##
## Marginal log-Likelihood: -10546.21
## 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 Relative risks
RR_stIntI.BYM <- c()
for (i in 1:316) {
RR_stIntI.BYM[i] <- inla.emarginal(
function(x) exp(x),
stIntI.FE.int.BYM.model$marginals.random$ID.space[[i]]
)
}
# Posterior probabilities (for spatial RR)
RR_stIntI.BYM_marg <- stIntI.FE.int.BYM.model$marginals.random$ID.space[1:316]
PP_stIntI.BYM <- lapply(RR_stIntI.BYM_marg, function(x) {
1 - inla.pmarginal(0, x)
})
# Temporal Relative risks and CI95
RR_stIntI.RW_RR <- c()
RR_stIntI.RW_lo <- c()
RR_stIntI.RW_hi <- c()
for (i in 1:7) {
# Posterior mean
RR_stIntI.RW_RR[i] <- inla.emarginal(
function(x) exp(x),
stIntI.FE.int.BYM.model$marginals.random$ID.time[[i]]
)
# 2.5% quantile
RR_stIntI.RW_lo[i] <- inla.qmarginal(0.025, inla.tmarginal(function(x) exp(x), stIntI.FE.int.BYM.model$marginals.random$ID.time[[i]]))
# 97.5% quantile
RR_stIntI.RW_hi[i] <- inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stIntI.FE.int.BYM.model$marginals.random$ID.time[[i]]))
}
RR_stIntI.RW <- data.frame(RR = RR_stIntI.RW_RR, low = RR_stIntI.RW_lo, high = RR_stIntI.RW_hi)
Temp2 <- ggplot(RR_stIntI.RW, aes(seq(2013, 2019), RR)) +
geom_line() +
ggtitle("ST model Int I") +
geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.2) +
labs(x = "year")
# print(Temp1)
print(Temp2)
ggsave("temp_RR.png", Temp2)
resRR_PP_stIntI <- data.frame(
resRR = RR_stIntI.BYM,
PP = unlist(PP_stIntI.BYM),
SP_ID = cases_agg[, 1]
)
# breakpoints
resRR_PP_stIntI$resRRcat <- cut(resRR_PP_stIntI$resRR, breaks = c(
min(resRR_PP_stIntI$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP_stIntI$resRR)
), include.lowest = T)
resRR_PP_stIntI$PPcat <- cut(resRR_PP_stIntI$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_ST.IntI <- left_join(map, resRR_PP_stIntI, by = c("lad19cd" = "Level"))
p5 <- ggplot() +
geom_sf(data = map_RR_ST.IntI) +
aes(fill = resRRcat) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = "RR")) +
ggtitle("a)") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
)
p6 <- ggplot() +
geom_sf(data = map_RR_ST.IntI) +
aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma",
name = "PP ST model Int I",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = "top",
reverse = T
)
) +
guides(fill = guide_legend(title = "Exceedence probability")) +
ggtitle("b)") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
)
print(p5)
ggsave("spatial_RR.png",p5)
print(p6)
ggsave("spatial_exceedence.png",p6)
# posterior probability
RR_stIntI_FE.BYM_marg <- stIntI.FE.int.BYM.model$marginals.random$ID.space.time[1:2212]
PP_stIntI_FE.BYM <- lapply(RR_stIntI_FE.BYM_marg, function(x) {
1 - inla.pmarginal(0, x)
})
resRR_stIntI_FE <- data.frame(
PP = unlist(PP_stIntI_FE.BYM),
SP_ID = cases_agg[, 1],
Year = rep(1:7, each = 316)
)
# join to map
map_RR_stIntI_FE <- left_join(map, resRR_stIntI_FE, by = c("lad19cd" = "Level"))
# breakpoints
map_RR_stIntI_FE$PPcat <- cut(resRR_stIntI_FE$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
# plot
p7<-ggplot() +
geom_sf(data = map_RR_stIntI_FE, aes(fill = PPcat)) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = NULL)) +
theme(
text = element_text(size = 20),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.9,0.2)
) +
facet_wrap(~Year, ncol = 4, labeller = labeller(Year = c("1" = "2013", "2" = "2014", "3" = "2015", "4" = "2016", "5" = "2017", "6" = "2018", "7" = "2019"))) +
labs("")
print(p7)
ggsave("int_map_PP.png", p7)
#c<-map_RR_stIntI_FE[map_RR_stIntI_FE$PP >0,][, c(3, 12, 13)]
DATA_ST_FE$intI <- stIntI.FE.int.BYM.model$summary.random$ID.space.time$mean
DATA_ST_FE$intI_cat <- cut(DATA_ST_FE$intI, breaks = c(
-1, -0.05,
-0.01, 0.01, 0.05, 1
), include.lowest = T)
p8 <- ggplot() +
geom_sf(data = DATA_ST_FE, aes(fill = intI_cat)) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = NULL)) +
theme(
text = element_text(size = 20),
axis.text.x = element_blank(),
axis.text.y = element_blank()
) +
facet_wrap(~Year, ncol = 3, labeller = labeller(ID.year = c("1" = "2013", "2" = "2014", "3" = "2015", "4" = "2016", "5" = "2017", "6" = "2018", "7" = "2019"))) +
labs("")
print(p8)
ggsave("int_map.png", p8)