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(kableExtra)
library(viridis)
library(RColorBrewer)
library(INLA)load("data.RData")
load("map.RData")Fixed effects model
formula_fixed <- O ~ PM25_scaled + Greenspace_log_scaled + IMD_scaled
fixed.model <- inla(formula_fixed, family = "poisson", data = data, 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.666, Running = 0.198, Post = 0.0261, Total = 0.89
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.520 0.002 -5.524 -5.520 -5.517 -5.520 0
## PM25_scaled -0.176 0.002 -0.179 -0.176 -0.172 -0.176 0
## Greenspace_log_scaled -0.027 0.002 -0.031 -0.027 -0.024 -0.027 0
## IMD_scaled 0.084 0.002 0.081 0.084 0.088 0.084 0
##
## Watanabe-Akaike information criterion (WAIC) ...: 43214.35
## Effective number of parameters .................: 97.15
##
## Marginal log-Likelihood: -26022.44
## 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 effects model
# aggregate data
cases_agg <- data %>%
group_by(Level) %>%
summarise(Observed = sum(O), Expected = sum(E), Population = sum(Population)) %>%
mutate(log.Population = log(Population)) %>%
dplyr::rename(O = Observed, E = Expected)
ID <- seq(1, 316)
formula_BYM2 <- O ~ 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.7, Running = 0.231, Post = 0.0145, Total = 4.95
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.559 0.006 -5.572 -5.559 -5.546 -5.559 0
##
## Random effects:
## Name Model
## ID BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 9.520 1.126 7.422 9.479 11.846 8.364
## Phi for ID 0.884 0.064 0.727 0.896 0.972 0.887
##
## Watanabe-Akaike information criterion (WAIC) ...: 3260.79
## Effective number of parameters .................: 158.93
##
## Marginal log-Likelihood: -1987.78
## 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)')RR_sBYM <- c()
for (i in 1:316) {
RR_sBYM[i] <- inla.emarginal(
function(x) exp(x),
sBYM.model$marginals.random$ID[[i]]
)
}
# Posterior probabilities
RR_sBYM_marg <- sBYM.model$marginals.random$ID[1:316]
PP_sBYM <- lapply(RR_sBYM_marg, function(x) {
1 - inla.pmarginal(0, x)
})
resRR_PP <- data.frame(
resRR = RR_sBYM,
PP = unlist(PP_sBYM),
SP_ID = cases_agg[, 1]
)
resRR_PP$resRRcat <- cut(resRR_PP$resRR, breaks = c(
min(resRR_PP$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP$resRR)
), include.lowest = T)
# breakpoints
resRR_PP$PPcat <- cut(resRR_PP$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_PP <- left_join(map, resRR_PP, by = c("lad19cd" = "Level"))
p1 <- ggplot() +
geom_sf(data = map_RR_PP) +
aes(fill = resRRcat) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = "RR")) +
ggtitle("RR Spatial model") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
)
p2 <- ggplot() +
geom_sf(data = map_RR_PP) +
aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma", name = "PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = "top",
reverse = T
)
) +
ggtitle("PP Spatial model") +
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(p1)
print(p2)
Spatio-temporal effects model
DATA_ST <- merge(map, data[c("Level", "Year", "O", "E", "log.Population")], by.x = "lad19cd", by.y = "Level")
# Create the ID for year (time)
DATA_ST$ID.time <- DATA_ST$Year - 2012
# Create the ID for space
DATA_ST$ID.space <- rep(seq(1, 316), each = 7)
formula_ST_noint <- O ~ 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)
)))
stBYM.model <- inla(
formula = formula_ST_noint, family = "poisson",
data = DATA_ST, offset = log.Population,
control.compute = list(waic = TRUE)
)
summary(stBYM.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.56, Running = 0.319, Post = 0.0167, Total = 4.9
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.564 0.007 -5.577 -5.564 -5.551 -5.564 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.517 1.126 7.412 9.478 11.835 8.357
## Phi for ID.space 0.885 0.064 0.729 0.897 0.972 0.887
## Precision for ID.time 205.903 132.125 47.663 174.469 545.847 7.694
##
## Watanabe-Akaike information criterion (WAIC) ...: 42104.54
## Effective number of parameters .................: 9342.73
##
## Marginal log-Likelihood: -13259.47
## 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_stBYM <- c()
for (i in 1:316) {
RR_stBYM[i] <- inla.emarginal(
function(x) exp(x),
stBYM.model$marginals.random$ID.space[[i]]
)
}
# Posterior probabilities (for spatial RR)
RR_stBYM_marg <- stBYM.model$marginals.random$ID.space[1:316]
PP_stBYM <- lapply(RR_stBYM_marg, function(x) {
1 - inla.pmarginal(0, x)
})
# Temporal Relative risks and CI95
RR_stRW_RR <- c()
RR_stRW_lo <- c()
RR_stRW_hi <- c()
for (i in 1:7) {
# Posterior mean
RR_stRW_RR[i] <- inla.emarginal(
function(x) exp(x),
stBYM.model$marginals.random$ID.time[[i]]
)
# 2.5% quantile
RR_stRW_lo[i] <- inla.qmarginal(0.025, inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
# 97.5% quantile
RR_stRW_hi[i] <- inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
}
RR_stRW <- data.frame(RR = RR_stRW_RR, low = RR_stRW_lo, high = RR_stRW_hi)
Temp1 <- ggplot(RR_stRW, aes(seq(2013, 2019), RR)) +
geom_line() +
ggtitle("ST model No Int") +
geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.2) +
labs(x = "year")
Temp1
resRR_PP_st <- data.frame(
resRR = RR_stBYM,
PP = unlist(PP_stBYM),
SP_ID = cases_agg[, 1]
)
# breakpoints
resRR_PP_st$resRRcat <- cut(resRR_PP_st$resRR, breaks = c(
min(resRR_PP_st$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP_st$resRR)
), include.lowest = T)
resRR_PP_st$PPcat <- cut(resRR_PP_st$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_ST <- left_join(map, resRR_PP_st, by = c("lad19cd" = "Level"))
p3 <- ggplot() +
geom_sf(data = map_RR_ST) +
aes(fill = resRRcat) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = "RR")) +
ggtitle("RR ST model") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
)
p4 <- ggplot() +
geom_sf(data = map_RR_ST) +
aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma",
name = "PP ST model",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = "top",
reverse = T
)
) +
ggtitle("PP ST model") +
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(p3)
print(p4)
Spatio-temporal effects with 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_intI <- O ~ 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.BYM.model <- inla(
formula = formula_ST_intI, family = "poisson", data = DATA_ST_FE, offset = log.Population,
control.compute = list(dic = TRUE, waic = TRUE)
)# Spatial Relative risks
RR_stIntI.BYM <- c()
for (i in 1:316) {
RR_stIntI.BYM[i] <- inla.emarginal(
function(x) exp(x),
stIntI.BYM.model$marginals.random$ID.space[[i]]
)
}
# Posterior probabilities (for spatial RR)
RR_stIntI.BYM_marg <- stIntI.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.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.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.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)
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("RR ST model Int I") +
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
)
) +
ggtitle("PP ST model Int I") +
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)
print(p6)
Mixed effects model
DATA_ST_FE$ID.space.time <- seq(1, dim(DATA_ST_FE)[1])
formula_ST_intI_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)
))) +
f(ID.space.time, model = "iid", hyper = list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)
)))
stIntI.FE.BYM.model <- inla(
formula = formula_ST_intI_FE, family = "poisson", data = DATA_ST_FE, offset = log.Population,
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(stIntI.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.48, Running = 0.807, Post = 0.0341, Total = 5.32
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.578 0.006 -5.589 -5.578 -5.566 -5.578 0
## PM25_scaled 0.034 0.016 0.002 0.033 0.066 0.033 0
## Greenspace_log_scaled -0.015 0.014 -0.043 -0.015 0.013 -0.015 0
## IMD_scaled 0.093 0.015 0.063 0.093 0.123 0.093 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.583 1.056 8.668 10.527 12.821
## Phi for ID.space 0.947 0.038 0.845 0.958 0.989
## Precision for ID.time 187.798 147.951 31.814 148.315 575.807
## Precision for ID.space.time 31.696 1.349 29.132 31.664 34.441
## mode
## Precision for ID.space 9.583
## Phi for ID.space 0.958
## Precision for ID.time 5.140
## Precision for ID.space.time 30.004
##
## Deviance Information Criterion (DIC) ...............: 18784.87
## Deviance Information Criterion (DIC, saturated) ....: 4077.90
## Effective number of parameters .....................: 1781.39
##
## Watanabe-Akaike information criterion (WAIC) ...: 18511.32
## Effective number of parameters .................: 1098.23
##
## Marginal log-Likelihood: -10536.91
## 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)')# Create table of results using kable# Spatial Relative risks
RR_stIntI_FE.BYM <- c()
for (i in 1:316) {
RR_stIntI_FE.BYM[i] <- inla.emarginal(
function(x) exp(x),
stIntI.FE.BYM.model$marginals.random$ID.space[[i]]
)
}
# Posterior probabilities (for spatial RR)
RR_stIntI_FE.BYM_marg <- stIntI.FE.BYM.model$marginals.random$ID.space[1:316]
PP_stIntI_FE.BYM <- lapply(RR_stIntI_FE.BYM_marg, function(x) {
1 - inla.pmarginal(0, x)
})
# Temporal Relative risks and CI95
RR_stIntI_FE.RW_RR <- c()
RR_stIntI_FE.RW_lo <- c()
RR_stIntI_FE.RW_hi <- c()
for (i in 1:7) {
# Posterior mean
RR_stIntI_FE.RW_RR[i] <- inla.emarginal(
function(x) exp(x),
stIntI.FE.BYM.model$marginals.random$ID.time[[i]]
)
# 2.5% quantile
RR_stIntI_FE.RW_lo[i] <- inla.qmarginal(0.025, inla.tmarginal(function(x) exp(x), stIntI.FE.BYM.model$marginals.random$ID.time[[i]]))
# 97.5% quantile
RR_stIntI_FE.RW_hi[i] <- inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stIntI.FE.BYM.model$marginals.random$ID.time[[i]]))
}
RR_stIntI_FE.RW <- data.frame(RR = RR_stIntI_FE.RW_RR, low = RR_stIntI_FE.RW_lo, high = RR_stIntI_FE.RW_hi)
Temp2 <- ggplot(RR_stIntI.RW, aes(seq(2013, 2019), RR)) +
geom_line() +
ggtitle("ST model Int I FE") +
geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.2) +
labs(x = "year")
print(Temp2)
resRR_PP_stIntI.FE <- data.frame(
resRR = RR_stIntI.BYM,
PP = unlist(PP_stIntI.BYM),
SP_ID = cases_agg[, 1]
)
# breakpoints
resRR_PP_stIntI.FE$resRRcat <- cut(resRR_PP_stIntI.FE$resRR, breaks = c(
min(resRR_PP_stIntI.FE$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP_stIntI.FE$resRR)
), include.lowest = T)
resRR_PP_stIntI.FE$PPcat <- cut(resRR_PP_stIntI.FE$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_ST.IntI.FE <- left_join(map, resRR_PP_stIntI.FE, by = c("lad19cd" = "Level"))
p7 <- ggplot() +
geom_sf(data = map_RR_ST.IntI.FE) +
aes(fill = resRRcat) +
theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill = guide_legend(title = "RR")) +
ggtitle("RR ST model Int I w/ FE") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
)
p8 <- ggplot() +
geom_sf(data = map_RR_ST.IntI.FE) +
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
)
) +
ggtitle("PP ST model Int I") +
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(p7)
print(p8)
DATA_ST_FE$intI <- stIntI.FE.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)
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("")
High-incidence locations
Spatial distribution of expected child hospital admission frequencies from 2013 to 2019.
# create bins for fill colour
DATA_ST_FE$E.cut <- cut(
DATA_ST_FE$E,
breaks = c(0, 50, 100, 200, 400, 600, 1000, max(DATA_ST_FE$E[!is.na(DATA_ST_FE$E)])),
)
# return rows where DATA_ST_FE$Rate.cut is NA
DATA_ST_FE[is.na(DATA_ST_FE$Rate.cut), ]## Simple feature collection with 0 features and 24 fields
## Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
## Projected CRS: WGS 84 / Pseudo-Mercator
## [1] lad19cd objectid_1 lad19nm
## [4] lad19nmw bng_e bng_n
## [7] long lat st_areasha
## [10] st_lengths GlobalID Year
## [13] O E log.Population
## [16] PM25_scaled Greenspace_log_scaled IMD_scaled
## [19] geometry ID.time ID.space
## [22] ID.space.time intI intI_cat
## [25] E.cut
## <0 rows> (or 0-length row.names)ggplot() +
geom_sf(data = DATA_ST_FE, aes(fill = E.cut)) +
theme_bw() +
scale_fill_brewer(name = "Expected child hospitalization", palette = "Reds", labels = c("<50", "(50, 100]", "(100, 200]", "(200, 400]", "(400, 600]", "(600, 1000]", "(1000, 1126]")) +
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("")