In njtierney/yahtsee: Yet Another Hierachical Time Series Extension and Expansion

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(yahtsee)
library(inlabru)
library(dplyr)

# add group id information
malaria_africa_ex <- malaria_africa_ts %>% 
  add_group_id(year) %>% 
  add_group_id(country) %>% 
  add_group_id(who_region) %>% 
  add_group_id(who_subregion) %>% 
  mutate(id = row_number()) %>% 
  mutate_if(is.character, as.factor) %>% 
  tidyr::drop_na() %>% 
  mutate(month_num = 12 * (year - min(year) + month))

# years <- rep(2000:2003, each =12)
# months <- rep(1:12, 4)
# 12 * (years - min(years)) + months

dplyr::n_distinct(malaria_africa_ex$month_num)
summary(malaria_africa_ex$month_num)

# some notes:
  # date/time has to be a factor, character, or numeric
  # random effect terms cannot be in the fixed effects
model <- INLA::inla(
  pr ~ lower_age + 
    f(
      date_num,
      group = .country_id,
      model = "ar1",
      constr = FALSE
    ),
  family = "gaussian",
  data = malaria_africa_ex,
  verbose = TRUE
  # options = list(
  #   control.compute = list(config = TRUE),
  #   control.predictor = list(compute = TRUE, link = 1)
  # )
)

model <- bru(
  pr ~ lower_age + Intercept + 
    r_country(
      month_num,
      group = .country_id,
      model = "ar1",
      constr = FALSE
    ) +
    r_subregion(
      month_num,
      model = "ar1",
      group = .who_subregion_id,
      constr = FALSE
      ) + 
    r_region(
      month_num,
      model = "ar1",
      group = .who_region_id,
      constr = FALSE
      ),
  family = "gaussian",
  data = malaria_africa_ex,
  options = list(
    control.compute = list(config = TRUE),
    control.predictor = list(compute = TRUE, link = 1)
)
)

model

new_test_data <- malaria_africa_ts %>% 
  mutate_if(is.character, as.factor) %>% 
  tidyr::drop_na() %>% 
  as_tibble() %>% 
  add_group_id(date) %>% 
  mutate(month_num = 12 * (year - min(year) + month)) %>%
  tsibble::as_tsibble(key = country,
                      index = .date_id)

m <- fit_hts(
  formula = pr ~ lower_age + upper_age +
    # inputs are  the levels of hierarchy, in order of decreasing size
    hts(who_region, who_subregion, country),
  .data = new_test_data,
  special_index = month_num,
  family = "gaussian"
)

library(tictoc)
tic()
identical(class(m), class(model))
toc()

## ****************************************
## inla observed
# previously, pred.inla1
model_inla_obs <- inla(
  treat_act ~ trtseek + yrid + f(ID, model = "iid"),
  family = "binomial",
  Ntrials = treat_total,
  data = as.data.frame(prop_data_act),
  # A boolean variable if the internal GMRF approximations be stored. 
  control.compute = list(config = TRUE),
  # compute is a boolean variable; should the marginals for the linear predictor be computed?
  # link is 
  control.predictor = list(compute = TRUE, link = 1)
)



## ## with dummy vars for regions + bym + s-p (linear effect of time for each area)
res.detrend.regions.bymsp <-
  inla(
    emplogitact.d ~ 0 + is_act + trtseek.d + popurban.d + anc1cov.d +
      oopfrac.d + hcaccess.d + eduallagsx.d + sbacov.d + haqi.d + 
      actyears20.d + logthepc.d + 
      who_subregion + year 
    f(
      country.idsp1, 
      year.idsp1,
      model = "iid"
    ) +
      f(
        year.idx,
        model = "ar1",
        group = country.id,
        constr = FALSE
      ) +
      f(
        year.idy,
        model = "ar1",
        group = whosubregion.id,
        constr = FALSE
      ) +
      f(
        year.idz,
        model = "ar1",
        group = whoregion.id,
        constr = FALSE
      ),
    data = data,
    quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975),
    control.inla = list(
      int.strategy =
        "eb"
    ),
    scale = emplogitactvar,
    control.compute = list(config = TRUE),
    control.predictor = list(compute = TRUE)
  )