Chapter 12 Exposure assessment-over space and time
In this chapter we have seen the many ways in which the time can be added to space in order to characterize random exposure fields. From this chapter, the reader will have gained an understanding of the following topics:
- Additional power that can be gained in an epidemiological study by combining the contrasts in the process over both time and space while characterizing the stochastic dependencies across both space and time for inferential analysis.
- Criteria that good approaches to spatio–temporal modelling should satisfy.
- General strategies for developing such approaches.
- Separability and non-separability in spatio–temporal models, and how these could be characterize using the Kronecker product of correlation matrices.
- Examples of the use of spatio–temporal models in modelling environmental exposures.
Example IDK: Spatio-temporal model for the UK ozone data
NOTE
The code for the implementation of DLMs in Stan and Nimble was developed by Paritosh Kumar Roy. This code and other more complex DLM structures are available through his github.
Nimble
library(dplyr)
library(sf)
library(ggplot2)
library(nimble, warn.conflicts = FALSE)
library(nleqslv)
source("functions/FFBS_functions_nimble.R")
load("data/ozone/uk_o3.rda")
load("data/ozone/uk_temp.rda")
load("data/ozone/uk_wind.rda")
load("data/ozone/uk_pollutant_coords.rda")
load("data/ozone/uk_pollutant_date.rda")
ls()## [1] "coords" "date"
## [3] "ddlm_moms_invariant" "ddlm_multi_obs"
## [5] "ddlm_uoms" "ddlm_uous"
## [7] "hyperParameterInvGamma" "nim_bsample"
## [9] "nim_chol2inv" "nim_diag"
## [11] "nim_ffbs_moms" "nim_ffbs_moms_invariant"
## [13] "nim_ffbs_uoms" "nim_ffbs_uous"
## [15] "nim_ffbs_woodbury" "nim_kf_moms"
## [17] "nim_postfit_uoms" "nim_postfit_uous"
## [19] "nim_woodbury" "nimSummary"
## [21] "o3" "rdlm_moms_invariant"
## [23] "rdlm_multi_obs" "rdlm_uoms"
## [25] "rdlm_uous" "sampler_ffbs_moms_invariant"
## [27] "sampler_ffbs_uoms" "sampler_ffbs_uous"
## [29] "sampler_ffbseta" "temp"
## [31] "wind"Example12_O3_Nimble <- nimbleCode({
sigma ~ T(dt(mu = 0, sigma = 1, df = 1), 0, Inf)
tau ~ T(dt(mu = 0, sigma = 1, df = 1), 0, Inf)
phi_inv ~ dgamma(shape = aPhi, rate = bPhi)
phi <- 1 / phi_inv
Vt[1:n, 1:n] <-
(sigma ^ 2) * exp(-obs_dist_mat[1:n, 1:n] / phi) + (tau ^ 2) * identityMatrix(d = n)
sqrt_Wt_diag[1] ~ T(dt(mu = 0, sigma = 1, df = 1), 0, Inf)
sqrt_Wt_diag[2] ~ T(dt(mu = 0, sigma = 1, df = 1), 0, Inf)
sqrt_Wt_diag[3] ~ T(dt(mu = 0, sigma = 1, df = 1), 0, Inf)
Wt[1:3, 1:3] <- nim_diag(x = sqrt_Wt_diag[1:3] ^ 2)
yt[1:n, 1:Tt] ~ ddlm_multi_obs(
Ft = Ft[1:p, 1:n, 1:Tt],
Vt = Vt[1:n, 1:n],
Gt = Gt[1:p, 1:p],
Wt = Wt[1:p, 1:p],
m0 = m0[1:p],
C0 = C0[1:p, 1:p]
)
# one can sample theta after fitting the model, using the posterior samples
# otherwise, sampling here is okay.
theta[1:p, 1:(Tt + 1)] <-
nim_ffbs_woodbury(
yt = yt[1:n, 1:Tt],
Ft = Ft[1:p, 1:n, 1:Tt],
Vt = Vt[1:n, 1:n],
Gt = Gt[1:p, 1:p],
Wt = Wt[1:p, 1:p],
m0 = m0[1:p],
C0 = C0[1:p, 1:p]
)
})# NOTE: I need to ask Paritosh about the ddlm_multi_obs function, it is missing
# in the functions file
# Model specification
n <- dim(o3)[1]
Tt <- dim(o3)[2]
p <- 3
Ft <- array(0, dim = c(p,n,Tt))
str(Ft)
Ft[1,1:n,1:Tt] <- 1
Ft[2,1:n,1:Tt] <- temp
Ft[3,1:n,1:Tt] <- wind
Gt <- diag(p)
Gt
# initials
m0 <- c(mean(o3),0,0)
C0 <- diag(x=c(1e2,1,1))
sqrt_Wt_diag <- c(0.1,0.06,0.06)
Wt <- diag(x = sqrt_Wt_diag^2)
tau <- 1
sigma <- 2
coords_sf <- st_as_sf(coords, coords = c("longitude", "latitude")) |>
st_set_crs(4326)
obs_dist_mat <- st_distance(coords_sf)
obs_dist_mat <- units::set_units(obs_dist_mat, km)
obs_dist_mat <- units::set_units(obs_dist_mat, NULL)
obs_max_dist <- max(obs_dist_mat)
obs_max_dist
obs_med_dist <- median(obs_dist_mat)
obs_med_dist
phi <- obs_med_dist / 6
phi
Vt <-
(sigma ^ 2) * exp(-obs_dist_mat / phi) + diag(x = tau ^ 2, nrow = n, ncol = n)
hyperPars <-
nleqslv(
c(2, 1),
hyperParameterInvGamma,
lower = 10,
upper = obs_med_dist,
prob = 0.98
)
prior_phi <-
1 / rgamma(n = 1000,
shape = hyperPars$x[1],
rate = hyperPars$x[2])
const_list <-
list(
n = n,
Tt = Tt,
p = p,
m0 = m0,
C0 = C0,
aPhi = hyperPars$x[1],
bPhi = hyperPars$x[2]
)
dat_list <-
list(
yt = o3,
Ft = Ft,
Gt = Gt,
obs_dist_mat = obs_dist_mat
)
init_list <-
list(
tau = 0.1,
sigma = 1,
sqrt_Wt_diag = rep(0.01, p),
phi_inv = 6 / obs_max_dist
)
Rmodel <-
nimbleModel(
Example12_O3_Nimble,
constants = const_list,
data = dat_list,
inits = init_list
)
Rmodel$calculate()
# -4862795
Rmodel$initializeInfo()
Cmodel <- compileNimble(Rmodel, showCompilerOutput = FALSE)
conf <-
configureMCMC(Rmodel,
monitors = c("tau", "sigma", "phi", "sqrt_Wt_diag", "theta"))
conf$removeSampler(target = "sqrt_Wt_diag[1]")
conf$addSampler(
target = "sqrt_Wt_diag[1]",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$removeSampler(target = "sqrt_Wt_diag[2]")
conf$addSampler(
target = "sqrt_Wt_diag[2]",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$removeSampler(target = "sqrt_Wt_diag[3]")
conf$addSampler(
target = "sqrt_Wt_diag[3]",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$removeSampler(target = "phi_inv")
conf$addSampler(
target = "phi_inv",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$removeSampler(target = "sigma")
conf$addSampler(
target = "sigma",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$removeSampler(target = "tau")
conf$addSampler(
target = "tau",
type = "RW",
control = list(
log = TRUE,
adaptive = TRUE,
adaptInterval = 100,
adaptFactorExponent = 0.8,
scale = 0.5
) ,
silent = TRUE
)
conf$printSamplers(byType = TRUE)
conf$printSamplers(executionOrder = TRUE)
Rmcmc <- buildMCMC(conf)
Cmcmc <-
compileNimble(
Rmcmc,
project = Cmodel,
resetFunctions = TRUE,
showCompilerOutput = FALSE
)
niter <- 10000
nburnin <- 0.5 * niter
nthin <- 1
start_time <- Sys.time()
post_samples <-
runMCMC(
Cmcmc,
niter = niter,
nburnin = nburnin,
thin = nthin,
nchains = 2,
samplesAsCodaMCMC = TRUE
)
end_time <- Sys.time()
run_time <- end_time - start_time
run_timelibrary(coda)
post_summary <- nimSummary(post_samples)
post_summary[c('tau','sigma','phi','sqrt_Wt_diag[1]','sqrt_Wt_diag[2]','sqrt_Wt_diag[3]'),]
library(tidyverse)
library(magrittr)
library(tidybayes)
tidy_post_samples <- post_samples %>% tidy_draws()
tidy_post_samples
tidy_post_samples %>%
select(.chain,.iteration,.draw,tau,sigma,phi,starts_with("sqrt_Wt_diag")) %>%
gather(vars,value,-.chain,-.iteration,-.draw) %>%
ggplot(aes(x=.iteration,y=value)) +
geom_path(aes(color = factor(.chain)), size = 0.25, show.legend = FALSE) +
facet_wrap(~vars, nrow = 2, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.background = element_blank())
ggsave(filename = "./figures/traceplot_integrated_model_ukO3Data55Sites.png", height = 6, width = 12)
tidy_post_samples %>%
select(.chain, .iteration, .draw,
starts_with("theta[1, 1]"),
starts_with("theta[1, 2]"),
starts_with("theta[1, 3]"),
starts_with("theta[1, 4]")) %>%
gather(vars,value,-.chain,-.iteration,-.draw) %>%
ggplot(aes(x=.iteration,y=value)) +
geom_path(aes(color = factor(.chain)), size = 0.25,
show.legend = FALSE) +
facet_wrap(~vars, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.background = element_blank())
tidy_post_samples %>%
select(.chain, .iteration, .draw,
starts_with("theta[2, 1]"),
starts_with("theta[2, 2]"),
starts_with("theta[2, 3]"),
starts_with("theta[2, 4]")) %>%
gather(vars,value,-.chain,-.iteration,-.draw) %>%
ggplot(aes(x=.iteration,y=value)) +
geom_path(aes(color = factor(.chain)), size = 0.25,
show.legend = FALSE) +
facet_wrap(~vars, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.background = element_blank())
tidy_post_samples %>%
select(.chain, .iteration, .draw,
starts_with("theta[3, 1]"),
starts_with("theta[3, 2]"),
starts_with("theta[3, 3]"),
starts_with("theta[3, 4]")) %>%
gather(vars,value,-.chain,-.iteration,-.draw) %>%
ggplot(aes(x=.iteration,y=value)) +
geom_path(aes(color = factor(.chain)), size = 0.25,
show.legend = FALSE) +
facet_wrap(~vars, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.background = element_blank())
post_sum_theta <- as.data.frame(post_summary) %>% rownames_to_column() %>%
filter(str_detect(rowname, "theta")) %>%
separate(rowname, into=c("x1","x2"), sep = ",") %>%
mutate(component = as.numeric(gsub(".*?([0-9]+).*", "\\1", x1))) %>%
mutate(time = as.numeric(gsub(".*?([0-9]+).*", "\\1", x2))) %>%
select(component,time,q2.5,q50,q97.5)
str(post_sum_theta)
ggplot(data = post_sum_theta, aes(x = time)) +
geom_ribbon(aes(ymin = q2.5,ymax = q97.5), fill = "lightgray", alpha = 0.7)+
geom_path(aes(y = q50), col = "blue", size = 0.4) +
facet_wrap(~component, nrow = 1, scales = "free",
labeller = label_bquote(theta[.(component)])) +
ylab("") +
xlab("Time") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.background = element_blank())
ggsave(filename = "./figures/states_integrated_model_ukO3Data55Sites.png", height = 3, width = 12)