## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  warning = FALSE,
  message = FALSE
)

options(rmarkdown.html_vignette.check_title = FALSE)

## ----Load required libraries, eval=FALSE--------------------------------------
# # Load necessary package dependencies
# library(dplyr)        # Data manipulation
# library(tidyr)        # Data tidying
# library(tidyselect)   # Helpers for selecting variables programmatically
# library(rlang)        # Tools for tidy evaluation and non-standard evaluation in tidyverse code
# library(ggplot2)      # Data visualization: creating plots and graphs
# library(cowplot)      # Combining and arranging multiple ggplot2 plots into a single figure
# library(grDevices)    # Base R graphics device functions (e.g., color palettes, saving plots)
# library(RColorBrewer) # Predefined color palettes for plots
# library(colorspace)   # Advanced color space manipulation and palettes
# library(sf)           # Handling spatial vector data (simple features)
# library(spdep)        # Spatial dependence and autocorrelation analysis
# library(sn)           # Skew-normal and skew-t distributions (for modelling skewed data)
# library(INLA)         # Integrated Nested Laplace Approximation for Bayesian models
# library(GHRexplore)   # Exploratory analysis of health data
# library(dlnm)         # Distributed linear and nonlinear models
# 
# # Load GHRmodel
# library(GHRmodel)

## ----Load required libraries invisible, echo=FALSE----------------------------
# Load necessary package dependencies
library(dplyr)        # Data manipulation
library(tidyr)        # Data tidying
library(tidyselect)   # Helpers for selecting variables programmatically
library(rlang)        # Tools for tidy evaluation and non-standard evaluation in tidyverse code
library(ggplot2)      # Data visualization: creating plots and graphs
library(cowplot)      # Combining and arranging multiple ggplot2 plots into a single figure
library(grDevices)    # Base R graphics device functions (e.g., color palettes, saving plots)
library(RColorBrewer) # Predefined color palettes for plots
library(colorspace)   # Advanced color space manipulation and palettes
library(sf)           # Handling spatial vector data (simple features)
library(spdep)        # Spatial dependence and autocorrelation analysis
library(sn)           # Skew-normal and skew-t distributions (for modeling skewed data)
library(GHRexplore)   # Exploratory analysis of health data

# Load GHRmodel
library(GHRmodel)     

## ----Process the data set-----------------------------------------------------
# Load data 
data("dengue_MS")
df <- dengue_MS

# Create ID variables
df <- df |>  
  # Filter out the year 2000 
  filter(year > 2000) |>  
  # Create numeric IDs for year, month and various spatial units.
  mutate(
    year_id = as.numeric(as.factor(year)),
    month_id = as.numeric(as.factor(month)),
    spat_id = as.numeric(as.factor(micro_code)), 
    spat_meso_id = as.numeric(as.factor(meso_code))
  )

## ----Load map_MS--------------------------------------------------------------
# Load map included in package
data("map_MS")

# Create adjacency Matrix
nb <- spdep::poly2nb(map_MS)
g <- spdep::nb2mat(nb, style = "B")

## ----Creating a onebasis for tmin---------------------------------------------
# Place 3 internal knots at the 25th, 50th, and 95th percentiles of tmin. 
tminknots = quantile(df$tmin,c(25,50,95)/100, na.rm=T)

# Construct a one-dimensional basis for modeling the nonlinear effect of minimum temperature (tmin).
ob_tmin <- onebasis_inla(
  covariate   = df$tmin,     # Numeric vector of the covariate values to transform.
  fun         = "bs",        # Basis function type: B-spline.
  knots       = tminknots,   # Internal knots for the spline
  degree      = 3,           # Polynomial degree (cubic spline)
  basis_name  = "tmin"       # Prefix for the generated basis matrix column names.
)

# Display the first few rows of the basis matrix to inspect the transformation.
head(ob_tmin)

## ----Create lagged matrices for tmin and pdsi---------------------------------
# Lagged covariate matrix for tmin 
lagged_tmin <- lag_cov(data = df,
                       name = "tmin",
                       time = "date",        
                       lag = c(1:2),         # 1 to 2-month lags
                       group = "micro_code", # Independent time-series for each microregion
                       add = FALSE)

# Lagged covariate matrix for pdsi
lagged_pdsi <- lag_cov(data = df,
                       name = "pdsi",
                       time = "date",        
                       lag = c(1:2),         # 1 to 2-month lags
                       group = "micro_code", # Independent time-series for each microregion
                       add = FALSE)


## ----Creating a crossbasis for tmin, eval=FALSE-------------------------------
# # 2-dimensional nonlinear effect of dengue risk across tmin values and lags
# cb_tmin <- crossbasis_inla(
#   covariate = lagged_tmin,              # Matrix of lagged tmin values
#   lag = c(1, 2),                        # Lags from 1 to 2 months
#   basis_name = "tmin",                  # Prefix for output column names
#   argvar  = list(fun = "bs"),           # Polynomial degree (cubic spline)
#   arglag  = list(fun = "ns",            # Natural cubic splines for lag-response
#                  knots=logknots(2,      # Maximum lag (here 2)
#                                 nk=1))  # Number of knots
# 
# 
# )
# 
# # Inspect the resulting matrix
# dim(cb_tmin)
# colnames(cb_tmin)

## ----load cb_tmin-------------------------------------------------------------
cb_tmin <- readRDS(system.file("examples", "cb_tmin.rds", package = "GHRmodel"))

## ----Creating a crossbasis for pdsi-------------------------------------------
# 2-dimensional nonlinear effect of dengue risk across pdsi values and lags
cb_pdsi <- crossbasis_inla(
  covariate = lagged_pdsi,
  lag = c(1, 2),
  basis_name = "pdsi",
  argvar  = list(fun = "bs"), # Default B-spline for exposure-response
  arglag  = list(fun = "bs")  # Default B-spline for lag-response
)

## ----Define monthly and yearly random effects priors--------------------------
# Define Gamma priors for the precision of temporal random effects
prior_t <- list(prec = list(prior = 'loggamma', param = c(0.01, 0.01))) 

## ----Define spatial random effect priors--------------------------------------
# Define penalized complexity (PC) priors for spatial random effects using BYM2
prior_sp <- list(
  prec = list(prior = 'pc.prec', param = c(0.5 / 0.31, 0.01)),  # Precision of spatial effect
  phi  = list(prior = 'pc', param = c(0.5, 2 / 3))              # Mixing parameter: structured vs unstructured
)


## ----Write DLNM INLA-compatible model formulas--------------------------------
# Write INLA-compatible model formulas including DLNM terms
formulas_dlnm <- write_inla_formulas(
  outcome = "dengue_cases",
  covariates = list(
    c("ob_tmin"),                 # Model 1: includes one-basis object (minimum temperature)
    c("cb_tmin"),                 # Model 2: includes cross-basis for temperature
    c("cb_tmin", "cb_pdsi")       # Model 3: includes cross-basis for temperature and drought index
    
  ),
  # First random effect: monthly seasonality
  re1 = list(
    id        = "month_id",      # Time index for monthly effect
    model     = "rw1",           # First-order random walk model
    cyclic    = TRUE,            # Enforce cyclicity over months (Dec → Jan wrap-around)
    hyper     = "prior_t",       # Hyperprior object for precision
    replicate = "spat_meso_id"   # Separate seasonal pattern per mesoregion
  ),
  # Second random effect: inter-annual trend
  re2 = list(
    id    = "year_id",           # Year index
    model = "iid",               # Independent and identically distributed
    hyper = "prior_t"            # Hyperprior for temporal smoothness
  ),
  # Third random effect: spatial autocorrelation
  re3 = list(
    id    = "spat_id",           # Spatial unit identifier
    model = "bym2",              # BYM2 spatial model (structured + unstructured)
    graph = "g",                 # Pre-computed spatial adjacency graph
    hyper = "prior_sp"           # PC priors for BYM2 model parameters
  ),
)

# Convert formulas list into a GHRformulas object for model fitting
formulas_dlnm_ghr <- as_GHRformulas(formulas_dlnm)

# Output is a GHRformulas object
class(formulas_dlnm_ghr)

## ----Fit DLNM INLA model, eval = FALSE----------------------------------------
# # Fitting a DLNM model using a previously defined formula and basis objects
# # Note: vcov = TRUE is required for DLNM interpretation and plotting
# model_dlnm <-fit_models(
#   formulas = formulas_dlnm_ghr,
#   name ="mod",
#   data = df,
#   family ="nbinomial",
#   offset = "population",
#   control_compute = list (config = TRUE, vcov = TRUE),
#   pb = TRUE)
# 
# class(model_dlnm)
# 

## ----fit_models model_user readrds, echo=FALSE--------------------------------
model_dlnm <- readRDS(system.file("examples", "model_dlnm.rds", package = "GHRmodel"))

## ----model_dlnm_gof, echo=FALSE-----------------------------------------------
# Extract goodness of fit metrics as a data.frame
model_dlnm_gof <- model_dlnm$mod_gof 
options(knitr.kable.NA = '')
knitr::kable(head(model_dlnm_gof[, c(1:3)]), 
             caption = "DLNM terms in `model_dlnm`")

## ----cpred_tmin---------------------------------------------------------------
# Generate dengue risk predictions for a range of tmin values. 
cpred_ob_tmin <- crosspred_inla(
  model = model_dlnm,   # GHRmodels object
  basis = ob_tmin,      # One-basis object for the `tmin` variable
  mod_id = "mod2",      # Identifier to select the correct model component
  cen = 19,             # Temperature centered at the mean of tmin 19°C
  from = 15,            # Predictions are generated for temperature values above 15
  to = 24               # Predictions are generated for temperature values below 24
)

# Output a GHRcrosspred object
class(cpred_ob_tmin)

## ----plot_coef_crosspred, fig.width = 6, fig.height = 4-----------------------
# Plot the predicted dengue risk values against the tmin values
plot_coef_crosspred(
  cpred_ob_tmin,           # The cross-predicted object from crosspred_inla()
  type = "slices",         # Plot type: shows effect as a continuous curve
  var = seq(15,24),        # Range of tmin values to display on the x-axis
  line_color = "red",      # Color of the fitted effect line
  line_size = 0.8,         # Thickness of the effect line
  ribbon_color = "red",    # Color of the credible interval (ribbon)
  ribbon_alpha = 0.3,      # Transparency level of the ribbon (0 = invisible)
  title = "Effect of minimum temperature on dengue risk",  # Plot title
  xlab = "Minimum temperature",   # Label for the x-axis
  ylab = "Effect"                 # Label for the y-axis 
)

## ----Evaluate coefficients for cb_tmin----------------------------------------
# Generate dengue risk predictions for a range of tmin values and lags. 
cpred_cb_tmin <- crosspred_inla(
  models = model_dlnm, # GHRmodels object
  basis = cb_tmin,     # Cross-basis from tmin
  mod_id = "mod3",     # Model identifier
  cen = 19,            # Centering value (mean minimum temperature 20°C)
  from = 17,           # Start predicting at 
  to= 24
)

# Output a GHRcrosspred object
class(cpred_cb_tmin)

## ----Plot coefficients for cb_tmin slices by lag, fig.width = 7, fig.height = 4----
# Plot the predicted dengue risk values against the tmin values and lags
plot_coef_crosspred(
  crosspred = cpred_cb_tmin,  # Crosspred object with model predictions
  type = "slices",            # Plot lag-specific slices of exposure-response curves
  exp = TRUE,                 # Exponentiate the coefficients,
  lag = 1:2,                  # Display results for lags 1 through 6
  line_color = "red",         # Red color for the lines representing effect estimates
  line_size = 0.8,            # Line thickness set to 0.8 for better visibility
  ribbon_color = "red",       # Red shading for credible interval ribbons
  ribbon_alpha = 0.3,         # Set ribbon transparency to 30%
  title = "Effect of minimum temperature on dengue relative risk by lag",  # Main plot title
  xlab = "Mean minimum temperature exposure",    # Label for the x-axis (exposure variable)
  ylab = "Relative Risk (RR)" # Label for the y-axis (effect estimate scale)
)

## ----Plot coefficients for cb_tmin slices by temp values,  fig.width = 7, fig.height = 5----
# Plot the predicted dengue risk values against the tmin values and lags
plot_coef_crosspred(
  crosspred = cpred_cb_tmin,  # Crosspred object with model predictions
  type = "slices",            # Plot temperature-specific slices of exposure-response curves
  exp = TRUE,                 # Exponentiate the coefficients (to relative risk scale)
  var = c(20:24),             # Display results for temperature 19°C to 24°C
  line_color = "red",         # Red color for the lines representing effect estimates
  line_size = 0.8,            # Line thickness set to 0.8 for better visibility
  ribbon_color = "red",       # Red shading for credible interval ribbons
  ribbon_alpha = 0.3,         # Set ribbon transparency to 30%
  title = "Effect of minimum temperatures 19°C to 24°C on dengue relative risk by lag",  # Main plot title
  xlab = "Lag",               # Label for the x-axis (exposure variable)
  ylab = "Relative Risk (RR)" # Label for the y-axis (effect estimate scale)
)

## ----Plot coefficients for cb_tmin heatmap, fig.width = 6, fig.height = 5-----
# Plot the predicted dengue risk values against the tmin values and lags
plot_coef_crosspred(
  crosspred = cpred_cb_tmin,  # Crosspred object containing model predictions
  type = "heatmap",           # Create a 2D heatmap of effects over exposure and lag values
  exp = TRUE,                 # Exponentiate effects to show relative risks
  palette = "-RdBu",          # Use the inverted "RdBu" color palette for the heatmap (blue for lower values)
  title = "Effect of minimum temperature on dengue relative risk by lag",  # Main plot title
  xlab = "Mean minimum temperature exposure",    # Label for the x-axis (exposure variable)
)


## ----Plot coefficients for cb_tmin heatmap smooth, fig.width = 6, fig.height = 5----
# Generate predictions using a higher-resolution range of temperature values
cpred_cb_tmin_smooth <- crosspred_inla(
  models = model_dlnm,    # GHRmodels object
  basis = cb_tmin,        # Cross-basis from tmin
  mod_id = "mod3",        # Model identifier
  cen = 19,               # Centering value (mean minimum temperature 19°C)
  at = seq(17,24, by=0.05), # Increase the number of temperature points for prediction
  from = 17, 
  to= 24
)

# Plot a more granular heatmap with increased number of temperature values for prediction
# and increased number of lag interpolation values
plot_coef_crosspred(
  crosspred = cpred_cb_tmin_smooth,   # Crosspred object containing model predictions
  type = "heatmap",           # Create a 2D heatmap of effects over exposure and lag
  exp = TRUE,                 # Exponentiate effects to show relative risks
  n_lag_smooth = 200,         # Smooth across lag with 200 interpolation points
  palette = "-RdBu",          # Use the inverted "RdBu" color palette
  title = "Effect of minimum temperature on dengue relative risk by lag",  # Main plot title
  xlab = "Mean minimum temperature exposure",    # Label for the x-axis (exposure variable)
)

## ----Plot coefficients for cb_tmin overall, fig.width = 6, fig.height = 4-----
# Plot overall effect of tmin on dengue across lags
plot_coef_crosspred(crosspred = cpred_cb_tmin,
                    type = "overall",
                    exp= TRUE, # Exponentiate effects to show relative risks
                    line_color = "Red", # Red color for the lines representing effect estimates
                    line_size = 0.8, # Line thickness set to 0.8 for better visibility
                    ribbon_color = "Red", # Red shading for credible interval ribbons
                    ribbon_alpha = 0.3, # Set ribbon transparency to 30%
                    xlab = "Mean minimum temperature exposure", # Label for the x-axis
                    ylab = "Effect" # Label for the y-axis (effect estimate scale)
)

