bam options to help with convergence

R/model-wis.R

@@ -44,7 +44,7 @@ m.formula <- wis ~
   s(Trend, bs = "re") +
   # Location
   s(location, bs = "re") +
-  # Week
+  # Week * location
   s(time, by = location, k = 40) +
   # Horizon
   s(Horizon, k = 3, by = Model, bs = "sz") +
@@ -61,16 +61,20 @@ m.fit <- function(outcomes, m.formula) {
       formula = m.formula,
       data = m.data |> filter(outcome_target == outcome),
       family = gaussian(link = "log"),
-      control = gam.control(trace = TRUE)
+      method = "fREML",
+      control = gam.control(trace = TRUE),
+      discrete = TRUE
     )
   })
 }
 # Fit
+cat("--------fitting univariate models")
 m.fits_uni_type <- m.fit(outcomes, m.formula_uni_type)
 m.fits_uni_tgt <- m.fit(outcomes, m.formula_uni_tgt)
 m.fits_uni_model <- m.fit(outcomes, m.formula_uni_model)
-m.fits_full <- m.fit(outcomes, m.formula)
-
+cat("--------fitting joint model")
+m.fits_joint <- m.fit(outcomes, m.formula)
+cat("finished fitting")
 # --- Output handling ---
 # Extract estimates for random effects
 random_effects_uni <- map_df(
@@ -79,14 +83,17 @@ random_effects_uni <- map_df(
   .id = "outcome_target") |>
   mutate(model = "Unadjusted")
 
-random_effects <- map_df(m.fits_full, extract_ranef,
+random_effects_joint <- map_df(m.fits_joint,
+                         extract_ranef,
                         .id = "outcome_target") |>
-  mutate(model = "Adjusted") |>
-  bind_rows(random_effects_uni)
+  mutate(model = "Adjusted")
 
-checks <- map(m.fits_full, k.check)
-formula <- m.fits[[1]]$formula
+random_effects <- random_effects_joint |>
+  bind_rows(random_effects_uni)
 
+# Extract model checks
+checks <- map(m.fits_joint, k.check)
+formula <- m.fits_joint[[1]]$formula
 results <- list(
   effects = random_effects,
   checks = checks,
@@ -95,7 +102,7 @@ results <- list(
 
 saveRDS(results, here("output", "results.rds"))
 
-iwalk(m.fits, \(x, target) {
+iwalk(m.fits_joint, \(x, target) {
   p <- appraise(x)
   ggsave(here("plots", paste0("check_", target, ".pdf")), p)
 })