gertvv / gemtc

GeMTC R package: model generation for network meta-analysis
GNU General Public License v3.0
43 stars 25 forks source link

Optional Burn in for the models - mtc.run #79

Open frostwix opened 1 month ago

frostwix commented 1 month ago

We want to make a proposition to update the mtc.run function as it is not supporting a burn-in phase currently. The burn-in phase is often required by national agencies for all mcmc runs, and even for conjugates. We make the mtc.run update as minimalistic as posible and backward compatible. That is why for backward compatiblity the NEW ARGUMENT n.burnin is added as the last one, so mtc.run will works even if somebody used it without argument names mtc.run(value, value) instead of mtc.run(arg1 = value, arg2 = value). Please give us your thoughts on the update.

Adaptation and burnin are not the same and currently only the former can be set in mtc.run. Please read more about it here https://stackoverflow.com/a/78807748/5442527 (example where adaptation is skipped as all models are conjugates) and here https://stackoverflow.com/a/38875637/5442527. Summing up it is important to recognize that burn-in and adaptation refer to two different things and are implemented independently in JAGS/rjags. The adaptation will be run only if at least one parameter in the model requires it, whereas burn-in is run always.

Our update fix the https://github.com/gertvv/gemtc/issues/78 issue.

We did not add any unit tests as the current tests structure will demand a lot of copy pasting. Please feel free to propose something.

Our update is backward compatible and results from mtc.run are the same. The default value for n.burnin is 0. Proof that the results from mtc.run are the same.

remotes::install_github("gertvv/gemtc/gemtc@master", force = TRUE)

library(gemtc)

set.seed(42)
model <- mtc.model(parkinson)
# By default, the model will have 4 chains - generate a seed for each
seeds <- sample.int(4, n = .Machine$integer.max)
# Apply JAGS RNG settings to each chain
model$inits <- mapply(c, model$inits, list(
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[1]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[2]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[3]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[4])), SIMPLIFY=FALSE)

results <- mtc.run(model, thin=1)

saveRDS(results, "mtc_run_before.RDS")

remotes::install_github("frostwix/gemtc/gemtc@master", force = TRUE)

library(gemtc)

set.seed(42)
model <- mtc.model(parkinson)
# By default, the model will have 4 chains - generate a seed for each
seeds <- sample.int(4, n = .Machine$integer.max)
# Apply JAGS RNG settings to each chain
model$inits <- mapply(c, model$inits, list(
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[1]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[2]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[3]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[4])), SIMPLIFY=FALSE)

results <- mtc.run(model, thin=1)

saveRDS(results, "mtc_run_after.RDS")

identical(readRDS("mtc_run_before.RDS")$samples, readRDS("mtc_run_after.RDS")$samples)
# TRUE
identical(readRDS("mtc_run_before.RDS")$deviance, readRDS("mtc_run_after.RDS")$deviance)
# TRUE