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. > library(epigrowthfit) > 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 ... > > > ## object ############################################################## > > o.1f <- fitted(o.1, class = TRUE, se = TRUE) > o.1f.e <- data.frame(top = gl(2L, 20L, labels = c("log(r)", "log(c0)")), + ts = gl(10L, 2L, 40L, labels = LETTERS[1:10]), + window = gl(20L, 1L, 40L, labels = sprintf("window_%02d", 1:20)), + value = as.double(o.1[["tmb_out"]][["env"]][[".__egf__."]][["adreport"]][["value"]]), + se = as.double(o.1[["tmb_out"]][["env"]][[".__egf__."]][["adreport"]][["sd"]])) > attr(o.1f.e, "se") <- TRUE > attr(o.1f.e, "ns") <- 20L > attr(o.1f.e, "nt") <- 2L > class(o.1f.e) <- c("fitted.egf", "data.frame") > stopifnot(identical(o.1f, o.1f.e)) > > > ## confint ############################################################# > > o.1fc <- confint(o.1f, level = 0.95, class = TRUE) > o.1fc.e <- o.1f[, c("top", "ts", "window", "value")] > o.1fc.e[["ci"]] <- + `dimnames<-`(epigrowthfit:::wald(o.1f[["value"]], o.1f[["se"]], + level = 0.95), + list(NULL, c("2.5 %", "97.5 %"))) > attr(o.1fc.e, "level") <- 0.95 > class(o.1fc.e) <- c("confint.egf", "data.frame") > stopifnot(identical(o.1fc, o.1fc.e)) > > proc.time() user system elapsed 4.70 0.25 4.95