## 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")
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.
# 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
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…
# 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.
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.
# 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')"
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",…
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")
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
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)
# 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.
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"
)