## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 5,
                      fig.height = 4.2, fig.align = "center", dpi = 100)

## -----------------------------------------------------------------------------
library(SDALGCP2)
library(sf)

data(sdalgcp_data)
head(sdalgcp_data)

# crude standardised incidence ratio (SIR): observed / expected-at-overall-rate
rate <- sum(sdalgcp_data$cases) / sum(sdalgcp_data$pop)
sdalgcp_data$SIR <- sdalgcp_data$cases / (sdalgcp_data$pop * rate)

## ----data-maps, fig.width = 7, fig.height = 3.4-------------------------------
plot(sdalgcp_data["SIR"], main = "Crude SIR (the data)")
plot(sdalgcp_data["x1"],  main = "Covariate x1")

## ----fit----------------------------------------------------------------------
set.seed(2024)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 4000, burnin = 1000, thin = 5,
                                         reanchor = 1))
summary(fit)

## ----risk-maps----------------------------------------------------------------
plot(fit, "relative_risk")   # relative risk exp(d'beta + S)
plot(fit, "adjusted_rr")     # covariate-adjusted relative risk exp(S)

## ----exceedance-maps----------------------------------------------------------
plot(fit, "adjusted_rr_se")                  # standard error of the adjusted RR
plot(fit, "exceedance", threshold = 1.5)     # P(adjusted RR > 1.5)

## ----continuous---------------------------------------------------------------
pc <- predict(fit, type = "continuous", sampler = "laplace", cellsize = 1)
plot(pc, "adjusted_rr", bound = sdalgcp_data)

## ----modelcheck---------------------------------------------------------------
chk <- model_check(fit)
chk$moran   # residual Moran's I and its permutation p-value