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 ... > > .S3method("all.equal", "confint.egf", + function(target, current, ignore = NULL, ...) { + if (!is.null(ignore)) + attributes(target)[ignore] <- attributes(current)[ignore] <- + NULL + NextMethod() + }) > > > ## object ############################################################## > > o.1c.w <- confint(o.1, A = NULL, method = "wald", class = TRUE, + random = TRUE) > o.1c.p <- confint(o.1, A = NULL, method = "profile", class = TRUE, + top = "log(r)", subset = quote(country == "A" & wave == 1)) Profile value: 1032.019 Profile value: 1032.051 Profile value: 1032.317 Profile value: 1032.555 Profile value: 1032.7 Profile value: 1032.863 Profile value: 1032.951 Profile value: 1033.042 Profile value: 1033.138 Profile value: 1033.238 Profile value: 1033.342 Profile value: 1033.45 Profile value: 1033.562 Profile value: 1033.677 Profile value: 1033.796 Profile value: 1033.918 Profile value: 1034.043 Profile value: 1032.019 Profile value: 1032.05 Profile value: 1032.28 Profile value: 1032.469 Profile value: 1032.58 Profile value: 1032.701 Profile value: 1032.831 Profile value: 1032.971 Profile value: 1033.119 Profile value: 1033.274 Profile value: 1033.437 Profile value: 1033.521 Profile value: 1033.607 Profile value: 1033.694 Profile value: 1033.782 Profile value: 1033.873 Profile value: 1033.964 > o.1c.u <- confint(o.1, A = NULL, method = "uniroot", class = TRUE, + top = "log(r)", subset = quote(country == "A" & wave == 1)) > > o.1f <- fitted(o.1, class = TRUE, se = TRUE) > o.1fc <- confint(o.1f, class = TRUE) > > o.1p <- profile(o.1, A = NULL, + top = "log(r)", subset = quote(country == "A" & wave == 1)) Profile value: 1032.019 Profile value: 1032.051 Profile value: 1032.317 Profile value: 1032.555 Profile value: 1032.7 Profile value: 1032.863 Profile value: 1032.951 Profile value: 1033.042 Profile value: 1033.138 Profile value: 1033.238 Profile value: 1033.342 Profile value: 1033.45 Profile value: 1033.562 Profile value: 1033.677 Profile value: 1033.796 Profile value: 1033.918 Profile value: 1034.043 Profile value: 1032.019 Profile value: 1032.05 Profile value: 1032.28 Profile value: 1032.469 Profile value: 1032.58 Profile value: 1032.701 Profile value: 1032.831 Profile value: 1032.971 Profile value: 1033.119 Profile value: 1033.274 Profile value: 1033.437 Profile value: 1033.521 Profile value: 1033.607 Profile value: 1033.694 Profile value: 1033.782 Profile value: 1033.873 Profile value: 1033.964 > o.1pc <- confint(o.1p, class = TRUE) > > stopifnot(exprs = { + is.list(o.1c.w) + identical(oldClass(o.1c.w), c("confint.egf", "data.frame")) + all.equal(o.1c.w, o.1fc) + + is.list(o.1c.p) + identical(oldClass(o.1c.p), c("confint.egf", "data.frame")) + all.equal(o.1c.p, o.1pc) + + is.list(o.1c.u) + identical(oldClass(o.1c.u), c("confint.egf", "data.frame")) + all.equal(o.1c.u, o.1c.p, tolerance = 5e-06) + }) > > > ## parallel ############################################################ > > f <- + function(method, cores) + confint(o.1, A = NULL, method = "uniroot", class = TRUE, + top = "log(r)", subset = quote(country == "A" & wave == 1), + parallel = egf_parallel(method = method, cores = cores)) > > windows <- .Platform[["OS.type"]] == "windows" > stopifnot(exprs = { + all.equal(o.1c.u, f("multicore", if (windows) 1L else 2L)) + all.equal(o.1c.u, f("snow", 2L)) + }) starting worker pid=38368 on localhost:11149 at 05:48:05.622 starting worker pid=169800 on localhost:11149 at 05:48:05.628 > > > ## plot ################################################################ > > op <- par(mar = c(4.5, 4, 2, 1), oma = c(0, 0, 0, 0)) > plot(o.1c.w) > par(op) > > proc.time() user system elapsed 41.62 0.76 52.89