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)
> library(methods)
> 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.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 
> stopifnot(exprs = {
+ 	is.list(o.1p)
+ 	identical(oldClass(o.1p), c("profile.egf", "profile"))
+ 	length(o.1p) == 1L
+ 	identical(names(o.1p), "log(r), A, window_01")
+ 	identical(dim(o.1p), c(1L, 1L))
+ 	identical(dimnames(o.1p), list("A, window_01", "log(r)"))
+ 
+ 	is.list(o.1p[[1L]])
+ 	identical(oldClass(o.1p[[1L]]), "data.frame")
+ 	length(o.1p[[1L]]) == 2L
+ 	identical(names(o.1p[[1L]]), c("z", "par.vals"))
+ 
+ 	is.double(z <- o.1p[[1L]][["z"]])
+ 	!is.matrix(z)
+ 	min(abs(z)) == 0
+ 	prod(sign(range(z))) == -1
+ 
+ 	is.double(par.vals <- o.1p[[1L]][["par.vals"]])
+ 	is.matrix(par.vals)
+ 	ncol(par.vals) == 1L # for now
+ 	!is.unsorted(par.vals, strictly = TRUE)
+ 	par.vals[which.min(abs(z))] == coef(o.1)[1L]
+ 
+ 	is.factor    (attr(o.1p, "top"   ))
+ 	is.factor    (attr(o.1p, "ts"    ))
+ 	is.factor    (attr(o.1p, "window"))
+ 	is.data.frame(attr(o.1p, "frame" ))
+ 	is           (attr(o.1p, "A"     ), "dgCMatrix")
+ 	is.double    (attr(o.1p, "par"   ))
+ 	identical    (attr(o.1p, "level" ), 0.95)
+ })
> 
> 
> ## confint #############################################################
> 
> o.1pc <- confint(o.1p, level = 0.95, class = TRUE)
> n <- length(o.1p)
> 
> stopifnot(exprs = {
+ 	is.list(o.1pc)
+ 	identical(oldClass(o.1pc), c("confint.egf", "data.frame"))
+ 	length(o.1pc) == 5L
+ 	identical(names(o.1pc), c("top", "ts", "window", "value", "ci"))
+ 
+ 	all(vapply(o.1pc[c("top", "ts", "window")], is.factor, FALSE))
+ 	all(vapply(o.1pc[c("value", "ci"        )], is.double, FALSE))
+ 
+ 	is.vector(o.1pc[["value"]])
+ 	is.matrix(o.1pc[["ci"]])
+ 	identical(dim(o.1pc[["ci"]]), c(1L, 2L))
+ 	identical(dimnames(o.1pc[["ci"]]), list(NULL, c("2.5 %", "97.5 %")))
+ 	all(o.1pc[["ci"]][, 1L] < o.1pc[["value"]])
+ 	all(o.1pc[["ci"]][, 2L] > o.1pc[["value"]])
+ })
> 
> 
> ## parallel ############################################################
> 
> f <-
+ function(method, cores)
+ 	profile(o.1, A = NULL,
+ 	        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.1p, f("multicore", if (windows) 1L else 2L))
+ 	all.equal(o.1p, f("snow", 2L))
+ })
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 
starting worker pid=39580 on localhost:11599 at 05:52:09.279
starting worker pid=47792 on localhost:11599 at 05:52:09.281
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 
> 
> 
> ## plot ################################################################
> 
> plot(o.1p, type = "z^2", bty = "u", las = 1)
> 
> proc.time()
   user  system elapsed 
  25.68    0.46   39.07