Load

## Load Harmonized Data ----
data <- read.csv("harmonized_data.csv")
head(data)
##   X year month       date    province dengue_cases population   tasmax
## 1 1 2023     1 2023-01-01  Xaisomboun            1     114000 23.89677
## 2 2 2023     1 2023-01-01      Attapu           10     166000 30.78710
## 3 3 2023     1 2023-01-01   Champasak            5     772000 30.78710
## 4 4 2023     1 2023-01-01      Xekong            3     134000 30.63548
## 5 5 2023     1 2023-01-01     Salavan            9     457000 29.26774
## 6 6 2023     1 2023-01-01 Savannakhet           15    1102000 28.32581
##      tasmin      tas prlr
## 1  7.970968 15.93387  8.8
## 2 20.403226 25.59516  0.8
## 3 18.535484 24.66129 32.8
## 4 16.403226 23.51935  5.0
## 5 17.719355 23.49355 19.2
## 6 15.461290 21.89355  0.7
## Load shapefile
laos<- st_read("gadm41_LAO1_cleaned.shp")
## Reading layer `gadm41_LAO1_cleaned' from data source 
##   `/home/akawieck/Documents/projects/dhis2.personal/ghrsuite_laos/gadm41_LAO1_cleaned.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 18 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 100.0868 ymin: 13.90968 xmax: 107.635 ymax: 22.5004
## Geodetic CRS:  WGS 84
ggplot() +
  geom_sf(data = laos)

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

GHRmodel description

Model dengue in Laos

We learned from the exploration of the data that precipitation, minimum temperature and average temperature might be the most relevant variables affecting dengue incidence. We can also observe that there may be a relationship between dengue outbreaks and increased average temperatures several months before.

Prepare data for modeling in INLA

# Assign numeric IDs to non-numeric variables for INLA modeling

data<- data %>%
  dplyr::mutate(dplyr::across(c("year", "month", "province"), 
                ~ as.numeric(as.factor(.)), 
                .names = "{.col}_id"), 
         dplyr::across(c("dengue_cases", "population", 
"tasmax", "tasmin", "tas", "prlr"), 
                ~ as.numeric(.))) %>%
  dplyr::select(-X) %>%
  # person time = 100000 person-month
  mutate(dengue_incidence=(dengue_cases / population) * 100000 ) # calculate incidence

Create lagged covariates

Here we create covariates lagged between 1-8 months for each observation, goruping by province.

# Lag covariates and attach to the original data
data <- ghrmodel::lag_cov(data = data,
                var = c("tas", "prlr"),
                time = "date",
                lag = c(8),
                group = "province",
                full = TRUE) # Merge = TRUE the matrix is merged to the data

dplyr::glimpse(data)
## Rows: 1,944
## Columns: 30
## $ year             <int> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015,…
## $ month            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5,…
## $ date             <chr> "2015-01-01", "2015-02-01", "2015-03-01", "2015-04-01…
## $ province         <chr> "Attapu", "Attapu", "Attapu", "Attapu", "Attapu", "At…
## $ dengue_cases     <dbl> 13, 39, 36, 21, 46, 68, 166, 171, 67, 44, 13, 8, 14, …
## $ population       <dbl> 166000, 166000, 166000, 166000, 166000, 166000, 16600…
## $ tasmax           <dbl> 31.98065, 33.42500, 36.29032, 37.68000, 36.58710, 33.…
## $ tasmin           <dbl> 17.55161, 20.41786, 24.76129, 25.89333, 26.89677, 25.…
## $ tas              <dbl> 24.78387, 26.94286, 30.55484, 31.81667, 31.76129, 29.…
## $ prlr             <dbl> 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, 293.6, 282.7, 148.9,…
## $ year_id          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,…
## $ month_id         <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5,…
## $ province_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ dengue_incidence <dbl> 7.831325, 23.493976, 21.686747, 12.650602, 27.710843,…
## $ tas.l1           <dbl> NA, 24.78387, 26.94286, 30.55484, 31.81667, 31.76129,…
## $ tas.l2           <dbl> NA, NA, 24.78387, 26.94286, 30.55484, 31.81667, 31.76…
## $ tas.l3           <dbl> NA, NA, NA, 24.78387, 26.94286, 30.55484, 31.81667, 3…
## $ tas.l4           <dbl> NA, NA, NA, NA, 24.78387, 26.94286, 30.55484, 31.8166…
## $ tas.l5           <dbl> NA, NA, NA, NA, NA, 24.78387, 26.94286, 30.55484, 31.…
## $ tas.l6           <dbl> NA, NA, NA, NA, NA, NA, 24.78387, 26.94286, 30.55484,…
## $ tas.l7           <dbl> NA, NA, NA, NA, NA, NA, NA, 24.78387, 26.94286, 30.55…
## $ tas.l8           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, 24.78387, 26.94286, 3…
## $ prlr.l1          <dbl> NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, 293.6, 282.7, 14…
## $ prlr.l2          <dbl> NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, 293.6, 282.7…
## $ prlr.l3          <dbl> NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, 293.6, 2…
## $ prlr.l4          <dbl> NA, NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, 293.…
## $ prlr.l5          <dbl> NA, NA, NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353.6, …
## $ prlr.l6          <dbl> NA, NA, NA, NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0, 353…
## $ prlr.l7          <dbl> NA, NA, NA, NA, NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 71.0,…
## $ prlr.l8          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, 7.6, 0.0, 2.3, 0.7, 7…

Explore non-linear effects

Priors

# Write  formulas with 3 random effects with customized prior
prior_re1 <- list(prec = list(prior = 'loggamma', param = c(0.01, 0.01)))
prior_re2 <- list(prec = list( prior = 'pc.prec', param = c(0.5 / 0.31, 0.01)),
                  phi = list( prior = 'pc', param = c(0.5, 2 / 3)))
Prior Type Parameters Meaning
prior_re1 Log-Gamma (0.01, 0.01) Weak prior on precision, allowing large variance.
prior_re2$prec PC Prior for Precision (0.5 / 0.31, 0.01) Shrinks precision toward a reasonable range, avoiding overfitting.
prior_re2$phi PC Prior for Spatial Dependency (0.5, 2/3) Encourages structured spatial correlation.

prior_re1: Log-Gamma Prior

The loggamma(0.01, 0.01) prior (shape = 0.01, rate = 0.01 ) is commonly used as a weakly informative prior, ensuring flexibility in the variance structure without forcing a strong assumption.

prior_re2: PC Prior for Precision

0.5 / 0.31 ≈ 1.61 sets the prior on precision.

0.01 is the probability that the standard deviation (σ = sqrt(1/precision)) exceeds a given threshold.

The Penalized Complexity (PC) prior is designed to avoid overfitting by shrinking unnecessary complexity. This implies a preference for a moderate variance but allows the data to override it if needed.

prior_re2: Spatial Dependency (phi)

0.5: The median of phi (controls spatial dependency balance).

2 / 3 ≈ 0.67: Probability that phi > 0.5.

phi is a mixing parameter that controls the balance between the structured spatial effect (neighbor-dependent) or unstructured spatial effect (independent random noise). A prior that favors phi > 0.5 encourages spatial smoothing, meaning that neighboring areas share more information.

Univariate predictors

First we will test the models with univariable predictors, to determine whether we should include linear or non-linear predictors in our model based on goodness of fit.

Write univariate model functions

# Create a list of linear lagged univariable predictors. 
# Includes all covariates that include the pattern tas.l and prlr.l 

cov_uni_l <- ghrmodel::extract_covariates(data=data, 
                                           pattern= c("tas.l", 
                                                      "prlr.l"))

dplyr::glimpse(cov_uni_l)
## List of 16
##  $ : chr "tas.l1"
##  $ : chr "tas.l2"
##  $ : chr "tas.l3"
##  $ : chr "tas.l4"
##  $ : chr "tas.l5"
##  $ : chr "tas.l6"
##  $ : chr "tas.l7"
##  $ : chr "tas.l8"
##  $ : chr "prlr.l1"
##  $ : chr "prlr.l2"
##  $ : chr "prlr.l3"
##  $ : chr "prlr.l4"
##  $ : chr "prlr.l5"
##  $ : chr "prlr.l6"
##  $ : chr "prlr.l7"
##  $ : chr "prlr.l8"
# Create a list of non linear lagged univariable predictors not replicated
cov_uni_l_nl <- ghrmodel::non_linear_covariates(covariates = cov_uni_l,
                                      method = "quantile",
                                      pattern = c("tas", "prlr"),
                                      n = 10)
dplyr::glimpse(cov_uni_l_nl)
## List of 16
##  $ : chr "f(INLA::inla.group(tas.l1, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l2, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l3, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l4, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l5, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l6, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l7, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l8, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l1, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l2, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l3, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l4, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l5, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l6, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l7, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l8, method='quantile', n=10), model='rw2')"
# Create a list of non linear lagged univariable predictors replicated by province
cov_uni_l_nl_rep <- ghrmodel::non_linear_covariates(covariates = cov_uni_l,
                                      method = "quantile",
                                      pattern = c("tas", "prlr"),
                                      n = 10,
                                      replicate = "province")

dplyr::glimpse(cov_uni_l_nl_rep)
## List of 16
##  $ : chr "f(INLA::inla.group(tas.l1, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l2, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l3, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l4, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l5, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l6, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l7, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(tas.l8, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l1, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l2, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l3, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l4, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l5, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l6, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l7, method='quantile', n=10), model='rw2', replicate =province)"
##  $ : chr "f(INLA::inla.group(prlr.l8, method='quantile', n=10), model='rw2', replicate =province)"
cov_uni_list <- c( "tas",
                     "prlr",
                      cov_uni_l, 
                      cov_uni_l_nl) #, cov_uni_l_nl_rep

dplyr::glimpse(cov_uni_list)
## List of 34
##  $ : chr "tas"
##  $ : chr "prlr"
##  $ : chr "tas.l1"
##  $ : chr "tas.l2"
##  $ : chr "tas.l3"
##  $ : chr "tas.l4"
##  $ : chr "tas.l5"
##  $ : chr "tas.l6"
##  $ : chr "tas.l7"
##  $ : chr "tas.l8"
##  $ : chr "prlr.l1"
##  $ : chr "prlr.l2"
##  $ : chr "prlr.l3"
##  $ : chr "prlr.l4"
##  $ : chr "prlr.l5"
##  $ : chr "prlr.l6"
##  $ : chr "prlr.l7"
##  $ : chr "prlr.l8"
##  $ : chr "f(INLA::inla.group(tas.l1, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l2, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l3, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l4, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l5, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l6, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l7, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(tas.l8, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l1, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l2, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l3, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l4, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l5, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l6, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l7, method='quantile', n=10), model='rw2')"
##  $ : chr "f(INLA::inla.group(prlr.l8, method='quantile', n=10), model='rw2')"

Write univariate model formulas

cov_uni_formulas <- write_inla_formulas(outcome = "dengue_cases",
                                        covariates = cov_uni_list ,
                                        re1 = list(id ="month_id",
                                                   re ="rw1", cyclic = TRUE,
                                                   hyper = "prior_re1",
                                                   replicate = "province_id" ),
                                        re2 = list(id = "year_id",
                                                   re = "rw1",
                                                   hyper = "prior_re1"),
                                        re3 = list(id = "province_id",
                                                   re = "bym2",
                                                   graph = "g", 
                                                   hyper = "prior_re2"),
                                                   baseline = TRUE)

head(cov_uni_formulas)
## [1] "dengue_cases ~ 1 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"         
## [2] "dengue_cases ~ 1 + tas + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"   
## [3] "dengue_cases ~ 1 + prlr + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"  
## [4] "dengue_cases ~ 1 + tas.l1 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"
## [5] "dengue_cases ~ 1 + tas.l2 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"
## [6] "dengue_cases ~ 1 + tas.l3 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1) + f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"
# transform formulas list into a GHRformulas object

cov_uni_formulas_ghr <- ghrmodel::as_GHRformulas(formulas = cov_uni_formulas)

class(cov_uni_formulas_ghr)
## [1] "GHRformulas" "list"
str(cov_uni_formulas_ghr)
## List of 4
##  $ formulas: chr [1:35] "dengue_cases ~ 1 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.mode"| __truncated__ "dengue_cases ~ 1 + tas + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scal"| __truncated__ "dengue_cases ~ 1 + prlr + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, sca"| __truncated__ "dengue_cases ~ 1 + tas.l1 + f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, s"| __truncated__ ...
##  $ vars    :'data.frame':    35 obs. of  1 variable:
##   ..$ covariate_1: chr [1:35] NA "tas" "prlr" "tas.l1" ...
##  $ re      : Named chr [1:3] "f(month_id, model = 'rw1', replicate = province_id, cyclic = TRUE, constr = TRUE, scale.model = TRUE, hyper = prior_re1)" "f(year_id, model = 'rw1', constr = TRUE, scale.model = TRUE, hyper = prior_re1)" "f(province_id, model = 'bym2', graph = g, constr = TRUE, scale.model = TRUE, hyper = prior_re2)"
##   ..- attr(*, "names")= chr [1:3] "re_1" "re_2" "re_3"
##  $ outcome : chr "dengue_cases"
##  - attr(*, "class")= chr [1:2] "GHRformulas" "list"
dplyr::glimpse(cov_uni_formulas_ghr$vars)
## Rows: 35
## Columns: 1
## $ covariate_1 <chr> NA, "tas", "prlr", "tas.l1", "tas.l2", "tas.l3", "tas.l4",…

Fit univariate models

m_uni <- ghrmodel::fit_models(formulas = cov_uni_formulas_ghr  ,
                           data = data,
                           family = "nbinomial",      # specify family
                           name = "m",
                           offset = "population",
                           config = TRUE,
                           pb = TRUE, 
                           nthreads = 8)
##   |                                           |                                   |   0%  |                                           |=                                  |   3% // Execution time: 0.15  min  |                                           |==                                 |   6% // Execution time: 0.25  min  |                                           |===                                |   9% // Execution time: 0.35  min  |                                           |====                               |  11% // Execution time: 0.46  min  |                                           |=====                              |  14% // Execution time: 0.56  min  |                                           |======                             |  17% // Execution time: 0.66  min  |                                           |=======                            |  20% // Execution time: 0.77  min  |                                           |========                           |  23% // Execution time: 0.88  min  |                                           |=========                          |  26% // Execution time: 0.98  min  |                                           |==========                         |  29% // Execution time: 1.08  min  |                                           |===========                        |  31% // Execution time: 1.19  min  |                                           |============                       |  34% // Execution time: 1.29  min  |                                           |=============                      |  37% // Execution time: 1.4  min  |                                           |==============                     |  40% // Execution time: 1.5  min  |                                           |===============                    |  43% // Execution time: 1.61  min  |                                           |================                   |  46% // Execution time: 1.72  min  |                                           |=================                  |  49% // Execution time: 1.83  min  |                                           |==================                 |  51% // Execution time: 1.94  min  |                                           |===================                |  54% // Execution time: 2.04  min  |                                           |====================               |  57% // Execution time: 2.18  min  |                                           |=====================              |  60% // Execution time: 2.3  min  |                                           |======================             |  63% // Execution time: 2.43  min  |                                           |=======================            |  66% // Execution time: 2.57  min  |                                           |========================           |  69% // Execution time: 2.71  min  |                                           |=========================          |  71% // Execution time: 2.84  min  |                                           |==========================         |  74% // Execution time: 2.96  min  |                                           |===========================        |  77% // Execution time: 3.09  min  |                                           |============================       |  80% // Execution time: 3.23  min  |                                           |=============================      |  83% // Execution time: 3.36  min  |                                           |==============================     |  86% // Execution time: 3.49  min  |                                           |===============================    |  89% // Execution time: 3.61  min  |                                           |================================   |  91% // Execution time: 3.74  min  |                                           |=================================  |  94% // Execution time: 3.86  min  |                                           |================================== |  97% // Execution time: 3.99  min  |                                           |===================================| 100% // Execution time: 4.12  min
# goodness of fit metrics in a dataframe
m_uni_gof <- m_uni$mod.gof

saveRDS(m_uni, "m_uni.rds")

Rank univariate models

WAIC

rank_uni_waic <- ghrmodel::rank_models(
  models = m_uni,
  metric = "waic",
  plot = TRUE,
  n = 8, 
  intercept = TRUE
)

rank_uni_waic_vs_base <- ghrmodel::rank_models(
  models = m_uni,
  metric = "waic_vs_base",
  plot = TRUE,
  n = 8, 
  intercept = TRUE
)

m_uni_gof %>%
  dplyr::filter(model_id %in% rank_uni_waic_vs_base)%>%
  dplyr::select(model_id, covariate_1, waic)%>%
  dplyr::arrange(waic)
##   model_id        covariate_1     waic
## 1      m33 prlr.l6_nl_q10_rw2 14187.03
## 2      m32 prlr.l5_nl_q10_rw2 14193.36
## 3      m35 prlr.l8_nl_q10_rw2 14194.57
## 4       m1               <NA> 14198.74
## 5      m23  tas.l4_nl_q10_rw2 14199.05
## 6      m34 prlr.l7_nl_q10_rw2 14199.07
## 7      m24  tas.l5_nl_q10_rw2 14199.80
## 8      m27  tas.l8_nl_q10_rw2 14200.50

CRPS

rank_uni_crps <- ghrmodel::rank_models(
  models = m_uni,
  metric = "crps",
  plot = TRUE,
  n = 8, 
  intercept = TRUE
)

m_uni_gof %>%
  dplyr::filter(model_id %in% rank_uni_crps)%>%
  dplyr::select(model_id, covariate_1, crps)%>%
  dplyr::arrange(crps)
##   model_id covariate_1     crps
## 1       m4      tas.l1 5.416200
## 2       m2         tas 5.425745
## 3       m8      tas.l5 5.443129
## 4       m5      tas.l2 5.498825
## 5      m10      tas.l7 5.515587
## 6       m9      tas.l6 5.523870
## 7      m11      tas.l8 5.531061
## 8       m6      tas.l3 5.539865

MAE

rank_uni_mae <- ghrmodel::rank_models(
  models = m_uni,
  metric = "mae",
  plot = TRUE,
  n = 8, 
  intercept = TRUE
)

m_uni_gof %>%
  dplyr::filter(model_id %in% rank_uni_mae)%>%
  dplyr::select(model_id, covariate_1, mae)%>%
  dplyr::arrange(mae)
##   model_id       covariate_1      mae
## 1      m20 tas.l1_nl_q10_rw2 50.13389
## 2      m22 tas.l3_nl_q10_rw2 50.47617
## 3      m21 tas.l2_nl_q10_rw2 50.56589
## 4      m11            tas.l8 50.96959
## 5      m23 tas.l4_nl_q10_rw2 51.05342
## 6       m8            tas.l5 51.08299
## 7       m7            tas.l4 51.09624
## 8       m9            tas.l6 51.12556

Posterior prediction check of univariate models

This function refits (or retrieves) the specified model, generates posterior predictions, and compares these predictions to the observed data by plotting their density estimates.

ppc_m1<-ghrmodel::post_pred_check(models = m_uni, model_id = "m1", s = 100, predictions = TRUE)

ppc_m33<-ghrmodel::post_pred_check(models = m_uni, model_id = "m33", s = 100, predictions = TRUE)

ppc_m4<-ghrmodel::post_pred_check(models = m_uni, model_id = "m4", s = 100, predictions = TRUE)

ppc_m20<-ghrmodel::post_pred_check(models = m_uni, model_id = "m20", s = 100, predictions = TRUE)

Fit vs. Observed plots for univariate models

# Plot model outputs
ghrmodel::plot_fit(
  models = m_uni,
  time = "date", 
  model_id = "m33",
  model_ref = "m1",
  #area = "province",
  selected_area = NULL,
  title = "Fitted (non-linear precipitation lag 6) vs Observed"
)

ghrmodel::plot_fit(
  models = m_uni,
  time = "date", 
  model_id = "m4",
  model_ref = "m1",
  #area = "province",
  selected_area = NULL,
  title = "Fitted (linear temperature lag 1) vs Observed"
)

ghrmodel::plot_fit(
  models = m_uni,
  time = "date", 
  model_id = "m20",
  model_ref = "m1",
  #area = "province",
  selected_area = NULL,
  title = "Fitted (non-linear temperature lag 1) vs. Observed"
)

Pretty similar to the only random effect model.

Plot random effects for univariate models

ghrmodel::plot_re_sp(
  models = m_uni,
  model_id = "m33",
  model_ref = "m1",
 # map = MS_map,
 #  map_area = "code",
  re_id = "province_id",
  label_model = NULL,
  title = "Spatial Random Effects"
)

ghrmodel::plot_re_sp(
  models = m_uni,
  model_id = "m33",
  model_ref = "m1",
  map = laos,
  map_area = "NAME_1",
  re_id = "province_id",
  label_model = NULL,
  title = "Spatial Random Effects"
)

ghrmodel::plot_re_t(
  models = m_uni,
  model_ids = c("m33", "m4", "m20"),
  model_ref = "m1",
  re_id = "year_id",
  label_model = NULL,
  title = "Yearly Random Effects"
)

ghrmodel::plot_re_t(
  models = m_uni,
  model_ids = c("m33", "m4", "m20"),
  model_ref = "m1",
  re_id = "month_id",
  replicated_id = "province_id",
  label_model = NULL,
  title = "Monthly Random Effects"
)