#' Old version of epilps
#'
#' @description
#' This routine estimates the instantaneous reproduction number Rt (the mean number
#' of secondary cases generated by an infectious individual at time t, White et
#' al. 2020) using Bayesian P-splines and Laplace approximations. Two
#' methods are possible for inference. LPSMAP is a fully sampling-free approach
#' based on Laplace approximations to the conditional posterior distribution of
#' the spline vector. LPSMALA is an MCMC-based approach based on Langevin
#' diffusions to sample the joint posterior of the model parameters. The
#' \code{epilps()} routine estimates Rt based on a time series of incidence
#' conts and a given serial interval distribution. The negative
#' Binomial distribution is used to model incidence count data and P-splines
#' (Eilers and Marx, 1996) are used to smooth the epidemic curve. The link
#' between the epidemic curve and the reproduction number is established via the
#' renewal equation. If \code{incidence} contains NA values at certain time
#' points, these are replaced by the average of the left- and right neighbor
#' counts. If the right neighbor is NA, the left neighbor is used as a
#' replacement value.
#'
#' @usage epilps(incidence, K = 30, method = c("LPSMAP","LPSMALA"),
#' serial_interval, penorder = 2, hyperprior = c(10,10),
#' chain_length = 5000, burn = 2000, progmala = TRUE, ci_level = 0.95,
#' etainit = c(1,5), cimethod = 1, verbose = TRUE, dates = NULL, tictoc = FALSE)
#'
#' @param incidence A vector containing the case counts per unit of time.
#' @param K Number of (cubic) B-splines in the basis.
#' @param method Either LPSMAP (fully sampling-free) or LPSMALA (MCMC-based).
#' @param serial_interval The discrete serial interval distribution.
#' @param penorder The order of the penalty (Default is second-order).
#' @param hyperprior Parameters for the Gamma prior on the dispersion parameter.
#' @param chain_length The length of the MCMC chain for method "LPSMALA"
#' (default 5,000).
#' @param burn The warm up period for method "LPSMALA" (default 2,000).
#' @param progmala Should the progress bar of LPSMALA be shown? (default TRUE).
#' @param ci_level Level of the credible intervals to be computed.
#' @param etainit Initial values for the hyperparameter vector (for the
#' optimization) in log scale.
#' @param cimethod The method used to construct credible intervals for Rt
#' with method LPSMAP. Default is 1 (log-normal approx) with scaling correction
#' on the covariance matrix. Setting it to 2 ignores the scaling correction.
#' @param verbose Should metainformation be printed?
#' @param dates A vector of date values (optional).
#' @param tictoc Should routine timing (in seconds) be measured?
#'
#' @return An object of class \code{epilps} containing the pointwise and set
#' estimates of the time-varying reproduction number and the epidemic curve
#' respectively.
#'
#' @author Oswaldo Gressani \email{oswaldo_gressani@hotmail.fr}
#'
#' @references Gressani, O., Wallinga, J., Althaus, C., Hens, N. and Faes, C.
#' (2021). EpiLPS: a fast and flexible Bayesian tool for near real-time
#' estimation of the time-varying reproduction number.
#' \emph{MedRxiv preprint}.
#' @references White, L.F., Moser, C.B., Thompson, R.N., Pagano, M. (2021).
#' Statistical estimation of the reproductive number from case notification
#' data. \emph{American Journal of Epidemiology}, \strong{190}(4):611-620.
#' @references Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing
#' with B-splines and penalties. \emph{Statistical Science},
#' \strong{11}(2):89-121.
#'
#' @examples
#' si <- c(0.05, 0.05, 0.1, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, 0.1, 0.1, 0.1)
#' epidemic <- episim(si = si, Rpattern = 2, endepi = 30)
#' # epifit <- epilps(incidence = epidemic$y, K = 30, serial_interval = si,)
#' # plot(epifit)
#'
#' @noRd
epilps <- function(incidence, K = 30, method = c("LPSMAP","LPSMALA"),
serial_interval, penorder = 2, hyperprior = c(10,10),
chain_length = 5000, burn = 2000, progmala = TRUE,
ci_level = 0.95, etainit = c(1,5), cimethod = 1, verbose = TRUE,
dates = NULL, tictoc = FALSE){
if (tictoc == TRUE) {
tic <- proc.time() # clock starts ticking
}
y <- incidence # time series of incidence counts (daily)
n <- length(y) # total number of days of the epidemic
smax <- length(serial_interval) # length of serial interval distribution
#-- Check incidence vector for NAs
if (anyNA(y)) {
nreplace <- sum(is.na(y))
NAprop <- round((nreplace / n) * 100, 2)
NAloc <- which(is.na(y))
for (j in 1:nreplace) {
if (1 < NAloc[j] && NAloc[j] < n) {
#NA is interior
if (!is.na(y[NAloc[j] + 1])) {
y[NAloc[j]] <- round((y[NAloc[j] - 1] + y[NAloc[j] + 1]) * 0.5)
} else{
y[NAloc[j]] <- y[NAloc[j] - 1]
}
} else if (NAloc[j] == 1) {
y[1] <- round(mean(y, na.rm = TRUE))
} else if (NAloc[j] == n) {
y[n] <- y[NAloc[j] - 1]
}
}
warning(paste0("Count data contains ", NAprop,
"% of NA values that have been replaced",
" (see documentation for details).",
" A 'too high' percentage of NA values can bias the analysis,",
" so care should be taken in interpreting the results in this case."))
}
#-- Assigning date vector
if (!is.null(dates)) {
datevec <- dates
} else{
datevec <- seq_len(n)
}
#-- B-splines basis
xx <- seq_len(n)
B <- Rcpp_KercubicBspline(xx, lower = 1, upper = n, K = K) # C++ call
#-- Penalty matrix
D <- diag(K)
for(k in 1:penorder){D <- diff(D)}
P <- t(D) %*% D # Penalty matrix of dimension c(K,K)
P <- P + diag(1e-06, K) # Perturbation to ensure P is full rank
#-- Parameter of Gamma prior for the roughness penalty parameter
a_delta <- hyperprior[1]
b_delta <- hyperprior[2]
nu <- 2
#-- Parameters of Gamma prior for the dispersion parameter
a_alpha <- 1e-04
b_alpha <- 1e-04
# log-likelihood, gradient and Hessian
loglik <- function(theta, a_disp) {
Btheta <- as.numeric(B %*% theta)
equal <- sum(lgamma(a_disp + y) - lgamma(a_disp)) +
sum(y * Btheta) + n * (a_disp * log(a_disp)) -
sum((y + a_disp) * log(a_disp + exp(Btheta)))
return(equal)
}
Dloglik <- function(theta, a_disp) {
Btheta <- as.numeric(B %*% theta)
res <- colSums((y - (y + a_disp) / (1 + a_disp * exp(-Btheta))) * B)
return(res)
}
D2loglik <- function(theta, a_disp) {
Btheta <- as.numeric(B %*% theta)
midvec <- a_disp * (y + a_disp) * (exp(Btheta) / (exp(Btheta) + a_disp) ^ 2)
Hess <- (-1) * (t(B) %*% (midvec * B))
return(Hess)
}
logptheta <- function(theta, a_disp, lambda) {
equal <- loglik(theta, a_disp) - 0.5 * lambda * sum((theta * P) %*% theta)
return(equal)
}
Dlogptheta <- function(theta, a_disp, lambda) {
equal <- Dloglik(theta, a_disp) - lambda * as.numeric(P %*% theta)
return(equal)
}
D2logptheta <- function(theta, a_disp, lambda) {
equal <- D2loglik(theta, a_disp) - lambda * P
return(equal)
}
# log-posterior of hyperparameter vector
logphyper <- function(vec_wv){
w <- vec_wv[1] # w = log(a_disp)
v <- vec_wv[2] # v = log(lambda)
# Alternative parameterization for Laplace routine
a_disp <- exp(w)
lambda <- exp(v)
# Laplace approximation
LL <- Rcpp_KerLaplace(theta0 = rep(1.5,K), a_disp, lambda, K, Dlogptheta, D2logptheta)
thetastar <- as.numeric(LL$Lapmode)
logdetSigstar <- Re(LL$logdetSigma)
equal <- 0.5 * logdetSigstar + 0.5 * (K + nu) * v + a_alpha * w -
(0.5 * nu + a_delta) * log(0.5 * nu * exp(v) + b_delta) +
loglik(thetastar, a_disp) -
0.5 * exp(v) * sum((thetastar * P) %*% thetastar) -
b_alpha * exp(w)
return(equal)
}
# optimization step
etamap <- stats::optim(etainit, fn=logphyper, control = list(fnscale=-1))$par
disphat <- exp(etamap[1])
lambhat <- exp(etamap[2])
Lap_approxx <- Rcpp_KerLaplace(theta0 = rep(1.5,K), disphat, lambhat, K, Dlogptheta, D2logptheta)
thetahat <- as.numeric(Lap_approxx$Lapmode)
Sighat <- Lap_approxx$Lapvar
#--- Metropolis-adjusted Langevin within Gibbs sampler
logtar <- function(zeta, lambda) {
theta <- zeta[1:K]
w <- zeta[(K + 1)]
rho <- exp(w)
equal <- loglik(theta, rho) -
0.5 * lambda * as.numeric(t(theta) %*% P %*% theta) -
b_alpha * exp(w) + a_alpha * w
return(equal)
}
Dlogtar <- function(zeta, lambda){
theta <- zeta[1:K]
w <- zeta[(K+1)]
rho <- exp(w)
Btheta <- as.numeric(B %*% theta)
grad_theta <- Dloglik(theta = theta, a_disp = rho) -
lambda * as.numeric(P %*% theta)
Dloglik_w <- sum(exp(w)*(digamma(y+exp(w))-digamma(exp(w))+(1+w)-
(log(exp(Btheta)+exp(w))+(1/(1+exp(Btheta-w)))))-
y * (1 / (1 + exp(Btheta - w))))
deriv_w <- Dloglik_w - b_alpha * exp(w) + a_alpha
res <- c(grad_theta, deriv_w)
return(res)
}
MALA <- function(M) {
SigLH <- matrix(0, nrow = (K + 1), ncol = (K + 1))
SigLH[1:K, 1:K] <- Sighat
SigLH[(K + 1), (K + 1)] <- 1
counter <- 0
zetamat <- matrix(0, nrow = M, ncol = (K + 1))
lambvec <- c()
deltavec <- c()
# Initial values
lambda_cur <- lambhat
w_cur <- etamap[1]
theta_cur <- thetahat
zeta_cur <- c(theta_cur, w_cur)
tun <- 0.15
#-- Progress bar
# if (progmala == TRUE) {
# progbar <- progress::progress_bar$new(
# format = crayon::cyan$bold("MALA running... [:bar] :percent"),
# total = M,
# clear = FALSE
# )
# }
for (m in 1:M) {
# New proposal
meanLH <- zeta_cur + 0.5 * tun * as.numeric(SigLH %*%
Dlogtar(zeta_cur, lambda_cur))
zeta_prop <- as.numeric(Rcpp_KerMVN(mu = meanLH, Sigma = (tun * SigLH)))
# Accept/Reject decision
G_cur <- Dlogtar(zeta_cur, lambda_cur)
G_prop <- Dlogtar(zeta_prop, lambda_cur)
ldiffq <- as.numeric((-0.5) * t(G_prop + G_cur) %*%
((zeta_prop - zeta_cur) +
(tun / 4) * SigLH %*% (G_prop - G_cur)))
ldiffp <- logtar(zeta_prop, lambda_cur) - logtar(zeta_cur, lambda_cur)
logr <- ldiffp + ldiffq
if (logr >= 0) {
zetamat[m,] <- zeta_prop
counter <- counter + 1
zeta_cur <- zeta_prop
} else if (logr < 0) {
u <- stats::runif(1)
if (u <= exp(logr)) {
zetamat[m,] <- zeta_prop
counter <- counter + 1
zeta_cur <- zeta_prop
} else{
zetamat[m, ] <- zeta_cur
}
}
# Gibbs step for delta
gdelta_shape <- 0.5 * nu + a_delta
gdelta_rate <- 0.5 * nu * lambda_cur + b_delta
deltavec[m] <- stats::rgamma(n = 1, shape = gdelta_shape,
rate = gdelta_rate)
# Gibbs step for lambda
glambda_shape <- 0.5 * (K + nu)
glambda_rate <- 0.5 * (as.numeric(t(zetamat[m,][1:K]) %*% P
%*% zetamat[m,][1:K]) + deltavec[m] * nu)
lambvec[m] <- stats::rgamma(n = 1, shape = glambda_shape,
rate = glambda_rate)
lambda_cur <- lambvec[m]
# Automatic tuning of Langevin algorithm
accept_prob <- min(c(1, exp(logr)))
heval <- sqrt(tun) + (1 / m) * (accept_prob - 0.57)
hfun <- function(x) {
epsil <- 1e-04
Apar <- 10 ^ 4
if (x < epsil) {
val <- epsil
} else if (x >= epsil && x <= Apar) {
val <- x
} else{
val <- Apar
}
return(val)
}
tun <- (hfun(heval)) ^ 2
# if (progmala == TRUE) {
# progbar$tick()
# }
}
accept_rate <- round(counter / M * 100, 2)
outlist <- list(theta_mcmc = zetamat[,(1:K)],
rho_mcmc = exp(zetamat[,(K+1)]),
lambda_mcmc = lambvec,
delta_mcmc = deltavec,
accept_rate = accept_rate,
niter = M)
return(outlist)
}
# Chosen method, LPSMAP or LPSMALA?
chosen_method <- match.arg(method)
if(chosen_method == "LPSMAP"){
# Point estimate and CI for mean of incidence counts
CImu <- function(t, alpha) {
bt <- as.numeric(Rcpp_KercubicBspline(t, lower = 1, upper = n, K = K))
zalpha <- stats::qnorm(1 - 0.5 * alpha, lower.tail = T)
sdd <- sqrt(sum((bt * Sighat) %*% bt))
logCIlow <- sum(thetahat * bt) - zalpha * sdd
logCIup <- sum(thetahat * bt) + zalpha * sdd
output <- cbind(exp(logCIlow), exp(logCIup))
return(output)
}
mu_estim <- as.numeric(exp(B %*% thetahat))
CI_mu <- t(sapply(seq_len(n), CImu, alpha = 1 - ci_level))
colnames(CI_mu) <- c(paste0("mu", ci_level * 100, "CI_low"),
paste0("mu", ci_level * 100, "CI_up"))
#----- Pointwise estimate of R(t) for LPS
Rt_LPS <- function(t) {
if (t == 1) {
res <- mu_estim[t]
} else if (t >= 2 && t <= smax) {
res <-
mu_estim[t] * ((sum(rev(mu_estim[1:(smax - 1)][1:(t - 1)]) *
serial_interval[1:(smax - 1)][1:(t - 1)])) ^ (-1))
} else if (t > smax && t <= n) {
res <- mu_estim[t] * (sum(rev(mu_estim[(t - smax):(t - 1)]) *
serial_interval) ^ (-1)) * (t > smax && t <= n)
}
return(res)
}
R_estim <- sapply(seq_len(n), Rt_LPS)
Rt_LPS_mean <- round(mean(R_estim[8:n]), 3)
muhat <- as.numeric(exp(B %*% thetahat))
Fano <- (1 + (1 / disphat) * mean(muhat)) / disphat
#------- (1-alpha) * 100 CI for R(t) at point t
CIRt_LPS <- function(t, alpha) {
hstar_t <- log(Rt_LPS(t))
muhat <- as.numeric(exp(B %*% thetahat))
dhstar_t <- c()
qz <- stats::qnorm(1 - 0.5 * alpha, lower.tail = TRUE)
for (k in 1:K) {
if (t == 1) {
add <- 0
} else if (t >= 2 && t <= smax) {
add <- (-1) * ((sum(rev(muhat[1:(smax - 1)][1:(t - 1)]) *
serial_interval[1:(smax - 1)][1:(t - 1)])) ^ (-1)) *
(sum(rev(muhat[1:(smax - 1)][1:(t - 1)]) *
serial_interval[1:(smax - 1)][1:(t - 1)] *
rev(B[, k][1:(smax - 1)][1:(t - 1)])))
} else if (t > smax && t <= n) {
add <- (-1) * (sum(rev(muhat[(t - smax):(t - 1)]) *
serial_interval) ^ (-1)) *
sum(rev(muhat[(t - smax):(t - 1)]) *
serial_interval * rev(B[, k][(t - smax):(t - 1)]))
}
dhstar_t[k] <- B[t, k] + add
}
if(cimethod == 1){
if (Fano <= 1) {
meanvar_ratio <- (1 / (1 + muhat[t] / disphat))
} else if (Fano > 1) {
meanvar_ratio <- 1
}
CIRt_low <- exp(hstar_t -
qz * sqrt(meanvar_ratio) *
sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t)))
CIRt_up <- exp(hstar_t +
qz * sqrt(meanvar_ratio) *
sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t)))
} else if (cimethod == 2){
CIRt_low <- exp(hstar_t -
qz * sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t)))
CIRt_up <- exp(hstar_t +
qz * sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t)))
}
# stats::qlnorm(alpha*0.5,meanlog=hstar_t,
# sdlog=((1/sqrt(1+muhat[t]/disphat)) *
# sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t))))
# stats::qlnorm(1-alpha*0.5,meanlog=hstar_t,
# sdlog=((1/sqrt(1+muhat[t]/disphat)) *
# sqrt(as.numeric(t(dhstar_t) %*% Sighat %*% dhstar_t))))
outlist <- list(CIRt_low = CIRt_low,
CIRt_up = CIRt_up)
return(outlist)
}
CI_R <- matrix(unlist(t(sapply(seq_len(n), CIRt_LPS, alpha = 1-ci_level))),
ncol = 2, byrow = FALSE)
colnames(CI_R) <- c(paste0("R", ci_level * 100, "CI_low"),
paste0("R", ci_level * 100, "CI_up"))
Rt_LPS_CImean <- round(colMeans(CI_R[8:n, ]), 3)
epifit <- data.frame(Date = datevec, R_estim, CI_R, mu_estim, CI_mu)
if (tictoc == TRUE) {
toc <- proc.time() - tic
toc <- round(toc[3], 3)
} else{
toc <- "Timer has not been requested."
}
outputlist <- list(CImu = CImu, Rt_LPS = Rt_LPS, CIRt_LPS = CIRt_LPS,
epifit = epifit, incidence = y,
serial_interval = serial_interval,
ci_level = ci_level, elapsed = toc, disphat = disphat)
}else if(chosen_method == "LPSMALA"){
Langevin <- MALA(M = chain_length)
# Chains and Geweke diagnostics
lambdaMALA <- coda::as.mcmc(Langevin$lambda_mcmc[(burn+1):chain_length])
deltaMALA <- coda::as.mcmc(Langevin$delta_mcmc[(burn+1):chain_length])
rhoMALA <- coda::as.mcmc(Langevin$rho_mcmc[(burn+1):chain_length])
thetaMALA <- coda::as.mcmc(Langevin$theta_mcmc[(burn+1):chain_length,])
Geweke <- matrix(0, nrow = (K + 3), ncol = 2)
rownames(Geweke) <- c(paste0("theta", seq_len(K)), "lambda", "delta", "rho")
colnames(Geweke) <- c("Geweke z-score", "Diagnostic passed (1=YES)")
Geweke[1:K, 1] <- unlist(coda::geweke.diag(thetaMALA))[1:K]
Geweke[(K + 1), 1] <- unlist(coda::geweke.diag(lambdaMALA))[1]
Geweke[(K + 2), 1] <- unlist(coda::geweke.diag(deltaMALA))[1]
Geweke[(K + 3), 1] <- unlist(coda::geweke.diag(rhoMALA))[1]
Geweke[, 2] <- as.logical(abs(Geweke[, 1]) < stats::qnorm(0.99))
# Estimated parameters
thetahat_mcmc <- colMeans(thetaMALA)
lambdahat_mcmc <- mean(lambdaMALA)
rhohat_mcmc <- mean(rhoMALA)
rhoCI <- coda::HPDinterval(rhoMALA, prob = ci_level)
accept_mcmc <- Langevin$accept_rate
muMALA_mcmc <- matrix(0, nrow = (Langevin$niter - burn), ncol = n)
for(j in 1:(Langevin$niter-burn)){
muMALA_mcmc[j,] <- as.numeric(exp(B %*% thetaMALA[j, ]))
}
mu_estim <- colMeans(muMALA_mcmc)
CI_mu <- coda::HPDinterval(coda::as.mcmc(muMALA_mcmc), prob = ci_level)
colnames(CI_mu) <- c(paste0("mu", ci_level * 100, "CI_low"),
paste0("mu", ci_level * 100, "CI_up"))
#----- Pointwise estimate of R(t) for MALA
Rt_MALA <- function(t) {
if (t == 1) {
res <- mu_estim[t]
} else if (t >= 2 && t <= smax) {
res <- mu_estim[t] *
((sum(rev(mu_estim[1:(smax - 1)][1:(t - 1)]) *
serial_interval[1:(smax - 1)][1:(t - 1)])) ^ (-1))
} else if (t > smax && t <= n) {
res <- mu_estim[t] * (sum(rev(mu_estim[(t - smax):(t - 1)]) *
serial_interval) ^ (-1)) * (t > smax && t <= n)
}
return(res)
}
R_estim <- sapply(seq_len(n), Rt_MALA)
Rt_LPSMALA_mean <- round(mean(R_estim[8:n]), 3)
#------- (1-alpha) * 100 CI for R(t) at point t with MALA
CIRt_MALA <- function(t, alpha) {
Rt_MALA_theta <- function(t, theta) {
fitted_mu_theta <- as.numeric(exp(B %*% theta))
if (t == 1) {
res <- fitted_mu_theta[t]
} else if (t >= 2 && t <= smax) {
res <- fitted_mu_theta[t] * ((sum(
rev(fitted_mu_theta[1:(smax - 1)][1:(t - 1)]) *
serial_interval[1:(smax - 1)][1:(t - 1)])) ^ (-1))
} else if (t > smax && t <= n) {
res <- fitted_mu_theta[t] *
(sum(rev(fitted_mu_theta[(t - smax):(t - 1)]) *
serial_interval) ^ (-1)) * (t > smax && t <= n)
}
return(res)
}
Rt_MALA_vec <- c()
for (j in 1:(Langevin$niter - burn)) {
Rt_MALA_vec[j] <- Rt_MALA_theta(t, thetaMALA[j, ])
}
Rt_MALA_vec <- coda::as.mcmc(Rt_MALA_vec)
HPDalpha <- coda::HPDinterval(Rt_MALA_vec, prob = (1 - alpha))
CIRt_low <- HPDalpha[1]
CIRt_up <- HPDalpha[2]
outlist <- list(CIRt_low = CIRt_low,
CIRt_up = CIRt_up)
return(outlist)
}
CI_R <- matrix(unlist(t(sapply(seq_len(n), CIRt_MALA,
alpha = 1 - ci_level))), ncol = 2, byrow = FALSE)
colnames(CI_R) <- c(paste0("R", ci_level * 100, "CI_low"),
paste0("R", ci_level * 100, "CI_up"))
Rt_LPSMALA_CImean <- round(colMeans(CI_R[8:n, ]), 3)
epifit <- data.frame(Date = datevec, R_estim, CI_R, mu_estim, CI_mu)
rownames(epifit) <- seq_len(n)
if (tictoc == TRUE) {
toc <- proc.time() - tic
toc <- round(toc[3], 3)
} else{
toc <- "Timer has not been requested."
}
outputlist <- list(chain_length = chain_length, burn = burn,
thetaMALA = thetaMALA, Rt_MALA = Rt_MALA,
CIRt_MALA = CIRt_MALA, epifit = epifit, K = K,
Geweke = Geweke, incidence = y,
serial_interval = serial_interval,
ci_level = ci_level, elapsed = toc,
disphat = rhohat_mcmc,
CIdisp = rhoCI)
}
if(verbose == TRUE){
if(chosen_method == "LPSMAP"){
cat("Inference method chosen: LPSMAP. \n", sep = "")
if(cimethod == 1){
cat("CI for LPSMAP computed via lognormal posterior approx. of Rt.")
} else if (cimethod == 2){
cat("CI for LPSMAP computed via sampling approach")
}
cat("Total number of days: ", n, ". \n",sep="")
cat("Mean Rt discarding first 7 days: ", Rt_LPS_mean, ".\n", sep="")
cat("Mean ",paste0(ci_level * 100,"%"),
" CI of Rt discarding first 7 days: (", Rt_LPS_CImean[1], ",",
Rt_LPS_CImean[2],") \n",sep="")
if (tictoc == TRUE) {
cat("Elapsed real time (wall clock time): ", toc, " seconds. \n", sep ="")
} else{
cat("Timing of routine not requested. \n")
}
} else if(chosen_method == "LPSMALA"){
cat("Inference method chosen: LPSMALA with chain length ", chain_length,
" and warmup ", burn, ".\n", sep="")
cat("MCMC acceptance rate: ", accept_mcmc, "%. \n", sep="")
cat("Geweke z-score < 2.33 for: ", sum(Geweke[,2]),"/",(K+3),
" variables. \n")
cat("Total number of days: ", n, ". \n",sep="")
cat("Mean Rt discarding first 7 days: ", Rt_LPSMALA_mean, ".\n", sep="")
cat("Mean ", paste0(ci_level * 100,"%"),
" CI of Rt discarding first 7 days: (", Rt_LPSMALA_CImean[1], ",",
Rt_LPSMALA_CImean[2],"). \n",sep="")
if (tictoc == TRUE) {
cat("Elapsed real time (wall clock time): ", toc, " seconds. \n", sep = "")
} else{
cat("Timing of routine not requested. \n")
}
}
}
attr(outputlist, "class") <- "epilps"
outputlist
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.