R Under development (unstable) (2024-06-17 r86768 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > attach(asNamespace("epigrowthfit")) > library(tools) > options(warn = 2L, error = if (interactive()) recover) > example("egf", package = "epigrowthfit"); o.1 <- m1; o.2 <- m2 egf> ## Simulate 'N' incidence time series exhibiting exponential growth egf> set.seed(180149L) egf> N <- 10L egf> f <- function(time, r, c0) { egf+ lambda <- diff(exp(log(c0) + r * time)) egf+ c(NA, rpois(lambda, lambda)) egf+ } egf> time <- seq.int(0, 40, 1) egf> r <- rlnorm(N, -3.2, 0.2) egf> c0 <- rlnorm(N, 6, 0.2) egf> data_ts <- egf+ data.frame(country = gl(N, length(time), labels = LETTERS[1:N]), egf+ time = rep.int(time, N), egf+ x = unlist(Map(f, time = list(time), r = r, c0 = c0))) egf> rm(f, time) egf> ## Define fitting windows (here, two per time series) egf> data_windows <- egf+ data.frame(country = gl(N, 1L, 2L * N, labels = LETTERS[1:N]), egf+ wave = gl(2L, 10L), egf+ start = c(sample(seq.int(0, 5, 1), N, TRUE), egf+ sample(seq.int(20, 25, 1), N, TRUE)), egf+ end = c(sample(seq.int(15, 20, 1), N, TRUE), egf+ sample(seq.int(35, 40, 1), N, TRUE))) egf> ## Estimate the generative model egf> m1 <- egf+ egf(model = egf_model(curve = "exponential", family = "pois"), egf+ formula_ts = cbind(time, x) ~ country, egf+ formula_windows = cbind(start, end) ~ country, egf+ formula_parameters = ~(1 | country:wave), egf+ data_ts = data_ts, egf+ data_windows = data_windows, egf+ se = TRUE) computing a Hessian matrix ... egf> ## Re-estimate the generative model with: egf> ## * Gaussian prior on beta[1L] egf> ## * LKJ prior on all random effect covariance matrices egf> ## (here there happens to be just one) egf> ## * initial value of 'theta' set explicitly egf> ## * theta[3L] fixed at initial value egf> m2 <- egf+ update(m1, egf+ formula_priors = list(beta[1L] ~ Normal(mu = -3, sigma = 1), egf+ Sigma ~ LKJ(eta = 2)), egf+ init = list(theta = c(log(0.5), log(0.5), 0)), egf+ map = list(theta = 3L)) computing a Hessian matrix ... > > > ## egf_has_random ###################################################### > > e <- new.env() > e[["data"]] <- list() > o <- list(tmb_out = list(env = e)) > class(o) <- "egf" > > e[["data"]][["Z"]] <- matrix(0, 3L, 3L) > stopifnot( egf_has_random(o)) > > e[["data"]][["Z"]] <- matrix(0, 3L, 0L) > stopifnot(!egf_has_random(o)) > > > ## egf_has_converged ################################################### > > o <- list(optimizer_out = list(convergence = 0L), + value = 100, + gradient = runif(10L, -0.5, 0.5), + hessian = TRUE) > class(o) <- c("egf", "list") # "list" only for dispatch to within.list > > stopifnot(exprs = { + egf_has_converged(o) + !egf_has_converged(within(o, optimizer_out[["convergence"]] <- 1L)) + !egf_has_converged(within(o, value <- NaN)) + !egf_has_converged(within(o, gradient[1L] <- 10)) + !egf_has_converged(within(o, hessian <- FALSE)) + is.na(egf_has_converged(within(o, hessian <- NA))) + }) > > > ## egf_par_expand #################################################### > ## egf_par_condense #################################################### > > t.2 <- o.2[["tmb_out"]] > p.2 <- t.2[["env"]][["last.par.best"]] > > p.2e <- egf_par_expand(t.2, p.2) > p.2ec <- egf_par_condense(t.2, p.2e) > > len <- c(table(factor(names(p.2), levels = c("beta", "theta", "b")))) > > stopifnot(exprs = { + all.equal(p.2e, structure(c(p.2[1:4], theta = 0, p.2[-(1:4)]), + len = len + c(0L, 1L, 0L), names = NULL)) + all.equal(p.2ec, structure(p.2, len = len, names = NULL)) + }) > > > ## egf_report ######################################################## > ## egf_adreport ######################################################## > > identical. <- + function(x, y, ...) { + x[["env"]] <- y[["env"]] <- NULL + identical(x, y, ...) + } > > r.1 <- o.1[["tmb_out"]][["env"]][[".__egf__."]][["report"]] > r.1. <- egf_report(o.1) > stopifnot(identical.(r.1., r.1)) > > o.1[["tmb_out"]][["env"]][[".__egf__."]][["report"]] <- NULL > r.1. <- egf_report(o.1) > stopifnot(identical.(r.1., r.1)) > > r.1 <- o.1[["tmb_out"]][["env"]][[".__egf__."]][["adreport"]] > r.1. <- egf_adreport(o.1) > stopifnot(identical.(r.1., r.1)) > > o.1[["tmb_out"]][["env"]][[".__egf__."]][["adreport"]] <- NULL > assertCondition(r.1. <- egf_adreport(o.1), "message") computing a Hessian matrix ... > stopifnot(identical.(r.1., r.1)) > > proc.time() user system elapsed 5.75 0.23 5.95