R/epilps.R

In EpiLPS: A Fast and Flexible Bayesian Tool for Estimating Epidemiological Parameters

#' 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

}