R/GDILM_SEIRS_Sim_Par_Est.R

In GDILM.SEIRS: Spatial Modeling of Infectious Disease with Reinfection

#' GDILM SEIRS for a Simulation Study
#'
#' This function conducts a simulation study for the Geographically Dependent Individual Level Model (GDILM) of infectious disease transmission, incorporating reinfection dynamics within the Susceptible-Exposed-Infectious-Recovered-Susceptible (SEIRS) framework, using a user-defined grid size. It applies a likelihood based Monte Carlo Expectation Conditional Maximization (MCECM) algorithm to estimate model parameters and compute the AIC.
#' @param GridDim1 First dimension of the grid
#' @param GridDim2 Second dimension of the grid
#' @param NPostPerGrid Number of postal codes per grid cell
#' @param MaxTimePand Last time point of the pandemic
#' @param tau0 Initial value for spatial precision
#' @param lambda0 Initial value for spatial dependence
#' @param alphaS0 Initial value for the susceptibility intercept
#' @param delta0 Initial value for the spatial decay parameter
#' @param alphaT0 Initial value for the infectivity intercept
#' @param PopMin Minimum population per postal code
#' @param PopMax Maximum population per postal code
#' @param InfFraction Fraction of each grid cell's population to be infected
#' @param ReInfFraction Fraction of each grid cell's population to be reinfected
#' @param InfPrd Infectious period that can be obtained either from the literature or by fitting an SEIRS model to the data
#' @param IncPrd Incubation period that can be obtained either from the literature or by fitting an SEIRS model to the data
#' @param NIterMC Number of MCMC iterations
#' @param NIterMCECM Number of MCECM iterations
#' @return
#'
#'   `alphaS` Estimate of alpha S
#'
#'   `BetaCovInf` Estimate of beta vector for the individual level infection covariate
#'
#'   `BetaCovSus` Estimate of beta vector for the areal susceptibility to first infection covariate
#'
#'   `BetaCovSusReInf` Estimate of beta vector for the areal susceptibility to reinfection covariate
#'
#'   `alphaT` Estimate of alpha T
#'
#'   `delta` Estimate of delta
#'
#'   `tau1` Estimate of tau
#'
#'   `lambda1` Estimate of lambda
#'
#'   `AIC` AIC of the fitted GDILM SEIRS
#'
#' @export
#' @import MASS
#' @import mvtnorm
#' @import ngspatial
#' @import stats
#'
#' @references Abed, A., Torabi, M., & Mashreghi, Z. (2025). Individual level modeling of infectious disease transmission with reinfection dynamics: Application to Tuberculosis in Manitoba, Canada. Spatial and Spatio-Temporal Epidemiology, 100780.
#'
#' @examples
#' \donttest{
#' # This example includes only one replication. The average of an arbitrary number
#' # of replications returns the estimation of the parameters.
#' # alphaS0 and alphaT0 must be small enough that the simulated epidemic does not
#' # infect every postal code immediately: with alphaS0 = 1 the whole grid is
#' # infected by t = 7, which leaves no susceptibles and no information to fit.
#' GDILM_SEIRS_Sim_Par_Est(5, 5, 10, 30, 0.7, 0.5, -12, 2.5, -3, 40, 50, 0.3, 0.6, 5, 5, 10, 3)
#' }
#'
GDILM_SEIRS_Sim_Par_Est=function(GridDim1,GridDim2,NPostPerGrid,MaxTimePand,tau0, lambda0, alphaS0, delta0, alphaT0,PopMin, PopMax,InfFraction,ReInfFraction, InfPrd, IncPrd, NIterMC, NIterMCECM){
  if(lambda0>=1|lambda0<0) stop("The spatial dependence parameter should be restricted to a range between 0 and 1.")
  if(lambda0==0) stop("Absence of spatial dependence: This model is designed for scenarios where spatial dependence is present.")
  if(delta0<=0) stop("The spatial decay parameter must be greater than zero.")
  if(NIterMC<=2) stop("The number of iterations must exceed 2.")
  if(InfPrd<=0) stop("The infectious period must be greater than zero.")
  if(IncPrd<=0) stop("The incubation period must be greater than zero.")
  if(InfFraction>1) stop("Fraction of each grid cell's population to be infected must be be restricted to a range between 0 and 1.")
  if(InfFraction<0) stop("Fraction of each grid cell's population to be infected must be be restricted to a range between 0 and 1.")
  if(ReInfFraction>1) stop("Fraction of each grid cell's population to be reinfected must be be restricted to a range between 0 and 1.")
  if(ReInfFraction<0) stop("Fraction of each grid cell's population to be reinfected must be be restricted to a range between 0 and 1.")
  NTotalGrid     <- GridDim1*GridDim2
  NAllPostPerGrid <- rep(NPostPerGrid, NTotalGrid)
  NTotalpost     <- sum(NAllPostPerGrid)
  # Place each cell's postal codes inside that cell's own square, walking the
  # cells in the same order `adjacency.matrix` indexes them. The previous version
  # blocked x by column and y by row against the same running index, which put a
  # whole run of cells in one square and left the geometry inconsistent with the
  # neighbourhood structure.
  generate_grid_data <- function(NPostPerGrid, GridDim1, GridDim2) {
    # adjacency.matrix(m, n) numbers its cells with the column varying fastest,
    # so walk the cells the same way to keep geometry and adjacency consistent.
    cells <- expand.grid(col = seq_len(GridDim2), row = seq_len(GridDim1))
    do.call(rbind, lapply(seq_len(nrow(cells)), function(g) {
      data.frame(xcor = runif(NPostPerGrid, cells$col[g]-1, cells$col[g]),
                 ycor = runif(NPostPerGrid, cells$row[g]-1, cells$row[g]))
    }))
  }

  data <- generate_grid_data(NPostPerGrid, GridDim1, GridDim2)
  Lat  <- data$xcor
  Long <- data$ycor
  AdjMat <- adjacency.matrix(GridDim1, GridDim2)
  D  <- (-1)*AdjMat
  diag(D) <- colSums(AdjMat)
  NLableGrid <- unlist(lapply(1:NTotalGrid, function(i) rep(i, NAllPostPerGrid[i])))
  NewLabelGrid <- sapply(1:NTotalGrid, function(g) rep(D[, g], NAllPostPerGrid))
  Dist <- as.matrix(dist(cbind(Lat, Long))) * 50
  Pop      <- sample(PopMin:PopMax, NTotalpost, replace = TRUE)
  NInf     <- pmax(round(Pop * InfFraction), 1)
  NReInf   <- pmax(round(Pop * ReInfFraction), 1)
  Sigma0 <- solve(tau0^2 * (lambda0 * D + (1-lambda0) * diag(NTotalGrid)))
  phi    <- mvrnorm(1, rep(0, NTotalGrid), Sigma0, tol = 1e-6)
  CovSus       <- cbind(rnorm(NTotalGrid,1,1), runif(NTotalGrid,1,2))
  CovSusReInf  <- cbind(rnorm(NTotalGrid,2,0.1), runif(NTotalGrid,2,3))
  CovInf       <- cbind(rnorm(NTotalpost,0,1), runif(NTotalpost,0,1))
  DimCovSus       <- ncol(CovSus)
  DimCovSusReInf  <- ncol(CovSusReInf)
  DimCovInf       <- ncol(CovInf)
  BetaCovSus0      <- rep(2, DimCovSus)
  BetaCovSusReInf0 <- rep(1, DimCovSusReInf)
  BetaCovInf0      <- rep(1, DimCovInf)
  ExpoTime  <- unlist(lapply(1:NTotalGrid, function(g) {
    vec <- rep(0, NAllPostPerGrid[g])
    vec[sapply(1:NTotalGrid, function(g) sample(NAllPostPerGrid[g], 1))[g]] <- 1
    vec
  }))

  ExpoTimeReInf <- unlist(lapply(1:NTotalGrid, function(g) {
    vec <- rep(0, NAllPostPerGrid[g])
    vec[sapply(1:NTotalGrid, function(g) sample(NAllPostPerGrid[g], 1))[g]] <- 1
    vec
  }))
  InfPeriod <- rep(InfPrd, NTotalpost)
  IncPeriod <- rep(IncPrd, NTotalpost)
  InfTime    <- ifelse(ExpoTime > 0, ExpoTime + IncPeriod[1], 0)
  ReInfTime  <- ifelse(ExpoTimeReInf > 0, ExpoTimeReInf + IncPeriod[1], 0)
  for(t in 1:MaxTimePand) {
    not_infected <- which(InfTime == 0)
    for(i in not_infected) {
      GridIndic <- NLableGrid[i]
      infectious <- which(NewLabelGrid[, GridIndic] != 0 &
                            InfTime <= t &
                            (InfTime + InfPeriod) >= t &
                            InfTime != 0)

      if(length(infectious) > 0) {
        dx <- sum(NInf[infectious] * exp(alphaT0 + CovInf[infectious, ] %*% BetaCovInf0) *
                    Dist[i, infectious]^(-delta0), na.rm = TRUE)
        P <- 1 - exp(-Pop[i] * exp(alphaS0 + CovSus[GridIndic, ] %*% BetaCovSus0 +
                                     CovSusReInf[GridIndic, ] %*% BetaCovSusReInf0 +
                                     phi[GridIndic]) * dx)
        if(runif(1) < P) InfTime[i] <- t + 1
      }
    }
  }
  # Keep the 0 sentinel for units that were never infected. Subtracting the
  # incubation period unconditionally made them negative, which satisfies neither
  # the susceptible branch nor is_exposed(), so they dropped out of the likelihood.
  ExpoTime <- ifelse(InfTime > 0, InfTime - IncPeriod[1], 0)
  for(t in 1:MaxTimePand) {
    not_infected <- which(ReInfTime == 0)
    for(i in not_infected) {
      GridIndic <- NLableGrid[i]
      infectious <- which(NewLabelGrid[, GridIndic] != 0 &
                            ReInfTime <= t &
                            (ReInfTime + InfPeriod) >= t &
                            ReInfTime != 0)
      if(length(infectious) > 0) {
        dx <- sum(NReInf[infectious] * exp(alphaT0 + CovInf[infectious, ] %*% BetaCovInf0) *
                    Dist[i, infectious]^(-delta0), na.rm = TRUE)
        P <- 1 - exp(-Pop[i] * exp(alphaS0 + CovSus[GridIndic, ] %*% BetaCovSus0 +
                                     CovSusReInf[GridIndic, ] %*% BetaCovSusReInf0 +
                                     phi[GridIndic]) * dx)
        if(runif(1) < P) ReInfTime[i] <- t + 1
      }
    }
  }
  ExpoTimeReInf <- ifelse(ReInfTime > 0, ReInfTime - IncPeriod[1], 0)
  # Fail fast, naming the parameter to change, instead of letting a rank-deficient
  # information matrix surface deep inside estfun() as a singular-matrix error.
  NInfected    <- sum(InfTime > 0)
  NSusceptible <- sum(ExpoTime == 0)
  NRegionsSusc <- length(unique(NLableGrid[ExpoTime == 0]))
  # One index case is seeded per grid cell, so anything at or below NTotalGrid means
  # the epidemic never actually transmitted.
  if(NInfected <= NTotalGrid){
    stop("The scenario setup doesn't include enough number of cases: the simulated",
         " epidemic produced ", NInfected, " infection(s) from ", NTotalGrid,
         " seeded index case(s), so there was no transmission to learn from. This",
         " function works only for a big number of cases; increase alphaS0, alphaT0,",
         " InfFraction or MaxTimePand.", call. = FALSE)
  }
  if(NSusceptible == 0){
    stop("The scenario setup doesn't include enough number of cases: the simulated",
         " epidemic infected all ", NTotalpost, " postal codes, so no susceptible unit",
         " is left and the area-level susceptibility information matrix is singular.",
         " This function works only for a big number of cases; use a smaller alphaS0",
         " and alphaT0 (for example alphaS0 = -12, alphaT0 = -3) so that the epidemic",
         " does not saturate.", call. = FALSE)
  }
  if(NRegionsSusc < DimCovSus){
    stop("The scenario setup doesn't include enough number of cases: susceptible units",
         " remain in only ", NRegionsSusc, " region(s), which cannot identify ",
         DimCovSus, " area-level susceptibility coefficients. This function works only",
         " for a big number of cases; use a larger grid, more postal codes per grid",
         " cell, or a smaller alphaS0.", call. = FALSE)
  }
  is_exposed <- function(ExpoTime, IncPeriod, t, i) {
    ExpoTime[i] <= t &
      (ExpoTime[i] + IncPeriod[i]) > t &
      ExpoTime[i] != 0
  }
  replace_nonfinite <- function(x, value = 0) {
    x[!is.finite(x)] <- value
    return(x)
  }
  # Stops when an information matrix carries too little information to be inverted
  # reliably, which happens when the simulated scenario contains too few cases.
  CheckEnoughCases <- function(M, what) {
    d <- det(M)
    if (!is.finite(d) || abs(d) < 1e-3) {
      stop("The scenario setup doesn't include enough number of cases: the ", what,
           " information matrix is numerically singular (|determinant| = ",
           format(abs(d), digits = 3),
           "). This function works only for a big number of cases; increase the number",
           " of infected and reinfected individuals and try again.", call. = FALSE)
    }
    invisible(d)
  }
  # Damped Newton step. Returns 0 (with a warning) instead of a meaningless update
  # when the curvature has collapsed numerically.
  NewtonStep <- function(score, hess, max_step = 1, tol = 1e-8) {
    if (!is.finite(score) || !is.finite(hess) || abs(hess) < tol) {
      warning("Degenerate Newton step: the curvature is numerically zero, so the",
              " parameter was left unchanged.", call. = FALSE)
      return(0)
    }
    # In every call site `hess` is `score` plus the second-order term. When that
    # term vanishes the model is saturated: every infection probability is
    # numerically 0 or 1, the step collapses to a constant and the parameter just
    # drifts by `max_step` per iteration instead of converging.
    if (abs(hess - score) <= tol*max(1, abs(score))) {
      warning("The model is saturated at the current parameter values: every infection",
              " probability is numerically 0 or 1, so the second-order term carries no",
              " information and the update reduces to a fixed step. Use a strongly",
              " negative alphaS0/alphaT0, or a scenario with more cases.", call. = FALSE)
    }
    max(min(score/hess, max_step), -max_step)
  }
  F1 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                   BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F11 <- rep(0, NTotalpost)
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
        F11[j] <- NInf[j] * exp(alphaT + CovInf[j,] %*% BetaCovInf) * Dist[i, j]^(-delta)
      }
    }
    F11<- replace_nonfinite(F11)
    return(sum(F11) * as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F2 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                        BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F22 <- rep(0, NTotalpost)
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && ReInfTime[j] <= t && (ReInfTime[j] + InfPeriod[j]) >= t && ReInfTime[j] != 0) {
        F22[j] <- NReInf[j] * exp(alphaT + CovInf[j,] %*% BetaCovInf) * Dist[i, j]^(-delta)
      }
    }
    F22<- replace_nonfinite(F22)
    return(sum(F22) * as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F3 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                   BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F33 <- array(0, c(DimCovInf, 1, NTotalpost))
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
        F33[,,j] <- NInf[j] * CovInf[j,] * as.numeric(exp(alphaT + CovInf[j,] %*% BetaCovInf)) * Dist[i, j]^(-delta)
      }
    }
    F33<- replace_nonfinite(F33)
    return(apply(F33, c(1,2), sum, na.rm = TRUE) *
             as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F4 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                       BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F44 <- array(0, c(DimCovInf, DimCovInf, NTotalpost))
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
        F44[,,j] <- NInf[j] * CovInf[j,] %*% t(CovInf[j,]) * as.numeric(exp(alphaT + CovInf[j,] %*% BetaCovInf)) * Dist[i, j]^(-delta)
      }
    }
    F44<- replace_nonfinite(F44)
    return(apply(F44, c(1,2), sum, na.rm = TRUE) *
             as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F5 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                       BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F55 <- rep(0, NTotalpost)
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
        F55[j] <- NInf[j] * exp(alphaT + CovInf[j,] %*% BetaCovInf) * Dist[i, j]^(-delta) * log(Dist[i, j])
      }
    }
    F55<- replace_nonfinite(F55)
    return(sum(F55) * as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F6 <- function(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                       BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F66 <- rep(0, NTotalpost)
    for (j in 1:NTotalpost) {
      if (NewLabelGrid[j, GridIndic] != 0 && InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
        F66[j] <- NInf[j] * exp(alphaT + CovInf[j,] %*% BetaCovInf) * Dist[i, j]^(-delta) * (log(Dist[i, j]))^2
      }
    }
    F66<- replace_nonfinite(F66)
    return(sum(F66) * as.numeric(exp(alphaS + CovSus[GridIndic,] %*% BetaCovSus + CovSusReInf[GridIndic,] %*% BetaCovSusReInf)))
  }
  F7 <- function(phi, alphaS, delta, lambda1, tau1, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    F8 <- array(0, c(NTotalpost, MaxTimePand, NTotalGrid))
    for (i in 1:NTotalpost) {
      for (t in 1:MaxTimePand) {
        for (GridIndic in 1:NTotalGrid) {
          if (NLableGrid[i] == GridIndic) {

            dx_fun <- function() {
              dx <- rep(0, NTotalpost)
              idx <- which(NewLabelGrid[, GridIndic] != 0 & InfTime <= t & (InfTime + InfPeriod) >= t & InfTime != 0)
              dx[idx] <- NInf[idx] * exp(alphaT + CovInf[idx, ] %*% BetaCovInf) * Dist[i, idx]^(-delta)
              sum(dx[is.finite(dx)])
            }
            if (ExpoTime[i] > t | ExpoTime[i] == 0) {
              dx <- dx_fun()

              F8[i, t, GridIndic] <- 1 - (1 - exp(-Pop[i] * exp(alphaS + CovSus[GridIndic, ] %*% BetaCovSus + phi[GridIndic] + CovSusReInf[GridIndic, ] %*% BetaCovSusReInf) * dx))
            }
            if (is_exposed(ExpoTime, IncPeriod, t, i)) {
              dx <- dx_fun()
              F8[i, t, GridIndic] <- 1 - exp(-Pop[i] * exp(alphaS + CovSus[GridIndic, ] %*% BetaCovSus + phi[GridIndic] + CovSusReInf[GridIndic, ] %*% BetaCovSusReInf) * dx)
            }
          }
        }
      }
    }
    # Drop the postal codes that contribute nothing (product of 1) with a
    # tolerance instead of round(., 10), which zeroed every product below 5e-11.
    F9 <- sapply(1:NTotalGrid, function(GridIndic) {
      F10 <- vapply(1:NTotalpost, function(i) {
        v <- F8[i, , GridIndic]
        prod(v[v > 0])
      }, numeric(1))
      F10 <- F10[is.finite(F10) & F10 > 0 & abs(F10 - 1) > 1e-10]
      if (length(F10) == 0) return(1)
      prod(F10)
    })
    return(F9)
  }
  alphaS=alphaS0
  delta=delta0
  tau1=tau0
  lambda1=lambda0
  BetaCovInf=BetaCovInf0
  BetaCovSus=BetaCovSus0
  BetaCovSusReInf=BetaCovSusReInf0
  alphaT=alphaT0
  Rnd <- matrix(0, NIterMC + 1, NTotalGrid)
  estfun=function(NLableGrid,Dist,alphaS,delta,lambda1,tau1,BetaCovInf,BetaCovSus,BetaCovSusReInf,alphaT){
    # Rebuild the proposal from the CURRENT variance components rather than from
    # the initial ones, so the sampler tracks the MCECM iterations.
    Sigma1 <- solve(tau1^2 * (lambda1 * D + (1 - lambda1) * diag(NTotalGrid)))
    Rnd[1, ] <- mvrnorm(1, rep(0, NTotalGrid), Sigma1, tol = 1e-6)
    for (L in 2:NIterMC) {
      phi <- mvrnorm(1, rep(0, NTotalGrid), Sigma1, tol = 1e-6)
      # Accept on the log scale: the product of per-region likelihoods underflows.
      LogNew <- sum(log(F7(phi, alphaS, delta, lambda1, tau1, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)))
      LogOld <- sum(log(F7(Rnd[L - 1, ], alphaS, delta, lambda1, tau1, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)))
      LogRatio <- LogNew - LogOld
      Rnd[L, ] <- if (is.finite(LogRatio) && log(runif(1)) < min(0, LogRatio)) phi else Rnd[L - 1, ]
    }
    Av1 <- function(Rnd, GridIndic) mean(exp(Rnd[1:NIterMC, GridIndic]))
    # Returns the escape probability q = 1 - P. `t` must be an argument: taking it
    # from the enclosing frame picked up a stale value left by the loops above.
    calc_esc <- function(L, i, GridIndic, t, Rnd, alphaS, delta, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
      idx <- which(NewLabelGrid[, GridIndic] != 0 & InfTime <= t & (InfTime + InfPeriod) >= t & InfTime != 0)
      contrib <- NInf[idx] * exp(alphaT + CovInf[idx, ] %*% BetaCovInf) * Dist[i, idx]^(-delta)
      dx <- sum(contrib[is.finite(contrib)])
      exp(-Pop[i] * exp(alphaS + CovSus[GridIndic, ] %*% BetaCovSus +
                          CovSusReInf[GridIndic, ] %*% BetaCovSusReInf + Rnd[L, GridIndic]) * dx)
    }
    Av2 <- function(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L) {
      q <- calc_esc(L, i, GridIndic, t, Rnd, alphaS, delta, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)
      if (!is.finite(q) || q >= 1) return(0)
      as.numeric(q / (1 - q) * exp(Rnd[L, GridIndic]))
    }
    Av3 <- function(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L) {
      q <- calc_esc(L, i, GridIndic, t, Rnd, alphaS, delta, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)
      if (!is.finite(q) || q >= 1) return(0)
      as.numeric(q / (1 - q)^2 * exp(2 * Rnd[L, GridIndic]))
    }
    S1 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        GridIndic <- NLableGrid[i]
        if (ExpoTime[i] > t | ExpoTime[i] == 0) {
          -Pop[i] * as.numeric(F1(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                                    BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
                                 Av1(Rnd, GridIndic))
        } else 0
      }), na.rm = TRUE)
    }), na.rm = TRUE)

    S2 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        GridIndic <- NLableGrid[i]
        if (is_exposed(ExpoTime, IncPeriod, t, i)) {
          SA4 <- sapply(1:NIterMC, function(L) {
            Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t,
                  BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)
          })
          Pop[i] * as.numeric(F1(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                                   BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
                                mean(SA4, na.rm = TRUE))
        } else 0
      }), na.rm = TRUE)
    }), na.rm = TRUE)
    S3 <- S1 + S2
    S4=rep(0,MaxTimePand)
    for(t in 1:MaxTimePand){
      S44=rep(0,NTotalpost)
      for(i in 1:NTotalpost){
        for(GridIndic in 1:NTotalGrid){
          if (NLableGrid[i]==GridIndic){
            if(is_exposed(ExpoTime, IncPeriod, t, i)){
              S444=c()
              for(L in 1:NIterMC){
                S444[L]=Av3(NLableGrid,Rnd,Dist,alphaS,delta,lambda1,i,GridIndic,t,BetaCovInf,BetaCovSus,BetaCovSusReInf,alphaT,L)
              }
              S44[i]=-Pop[i]^2*as.numeric((F1(NLableGrid,Dist,alphaS,delta,i,GridIndic,t,BetaCovInf,BetaCovSus,BetaCovSusReInf,alphaT))^2*mean(S444))
            }
          }
        }
      }
      S4[t]=sum(S44)
    }
    S5=sum(S4,na.rm=T)
    S6=S3+S5
    EstAlphaS=alphaS-NewtonStep(S3,S6)
    S7 <- array(0, c(DimCovSus, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(ExpoTime[i] > t | ExpoTime[i] == 0){
          S7[,,t] <- S7[,,t] - Pop[i] * CovSus[GridIndic, , drop = FALSE] *
            as.numeric(F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                            BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
                         Av1(Rnd, GridIndic))
        }
      }
    }
    S8 <- apply(S7, c(1,2), sum)
    S9 <- array(0, c(DimCovSus, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          S99 <- sapply(1:NIterMC, function(L)
            Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                  i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L))
          S9[,,t] <- S9[,,t] + Pop[i] * CovSus[GridIndic, , drop = FALSE] *
            as.numeric(F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                            BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
                         mean(S99))
        }
      }
    }
    S10 <- apply(S9, c(1,2), sum, na.rm = TRUE)
    S11 <- S8 + S10
    S12 <- array(0, c(DimCovSus, DimCovSus, MaxTimePand))
    S13 <- array(0, c(DimCovSus, DimCovSus, MaxTimePand))
    S14  <- array(0, c(DimCovSus, DimCovSus, MaxTimePand))
    CovMatrices <- lapply(1:NTotalGrid, function(g) {
      cov_vec <- matrix(CovSus[g, ], ncol = 1)
      cov_vec %*% t(cov_vec)
    })
    for (t in 1:MaxTimePand) {
      for (i in 1:NTotalpost) {
        GridIndic <- NLableGrid[i]
        if (ExpoTime[i] > t | ExpoTime[i] == 0) {
          contrib <- -Pop[i] * as.numeric(
            F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                 BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
              Av1(Rnd, GridIndic)
          )
          S12[,,t] <- S12[,,t] + CovMatrices[[GridIndic]] * contrib
        }
      }
      for (i in 1:NTotalpost) {
        GridIndic <- NLableGrid[i]
        if (is_exposed(ExpoTime, IncPeriod, t, i)) {
          S15 <- numeric(NIterMC)
          for (L in 1:NIterMC) {
            S15[L] <- Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                             i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)
          }
          contrib <- Pop[i] * as.numeric(
            F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                 BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) *
              mean(S15)
          )
          S13[,,t] <- S13[,,t] + CovMatrices[[GridIndic]] * contrib
        }
      }
      for (i in 1:NTotalpost) {
        GridIndic <- NLableGrid[i]
        if (is_exposed(ExpoTime, IncPeriod, t, i)) {
          S16 <- numeric(NIterMC)
          for (L in 1:NIterMC) {
            S16[L] <- Av3(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                             i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)
          }
          contrib <- -Pop[i]^2 * as.numeric(
            (F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                  BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)^2) * mean(S16)
          )
          S14[,,t] <- S14[,,t] + CovMatrices[[GridIndic]] * contrib
        }
      }
    }
    S17 <- apply(S12, c(1,2), sum)
    S18 <- apply(S13, c(1,2), sum, na.rm = TRUE)
    S19  <- apply(S14, c(1,2), sum, na.rm = TRUE)
    S20  <- S17 + S18 + S19 + 0.01
    CheckEnoughCases(S20, "area-level susceptibility")
    EstBetaCovSus <- BetaCovSus - solve(S20) %*% S11
    S21 <- array(0, c(DimCovSusReInf, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      S22 <- array(0, c(DimCovSusReInf, 1, NTotalpost))
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(ExpoTime[i] > t | ExpoTime[i] == 0){
          contrib <- -Pop[i] * CovSusReInf[GridIndic,] *
            as.numeric(F2(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                 BetaCovInf, EstBetaCovSus, BetaCovSusReInf, alphaT) *
                         Av1(Rnd, GridIndic))
          S22[,,i] <- contrib
        }
      }
      S21[,,t] <- apply(S22, c(1,2), sum)
    }
    S23 <- apply(S21, c(1,2), sum, na.rm = TRUE)
    S24 <- array(0, c(DimCovSusReInf, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      S25 <- array(0, c(DimCovSusReInf, 1, NTotalpost))
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          S26 <- sapply(1:NIterMC, function(L) Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                                                      i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L))
          contrib <- Pop[i] * CovSusReInf[GridIndic,] *
            as.numeric(F2(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                 BetaCovInf, EstBetaCovSus, BetaCovSusReInf, alphaT) *
                         mean(S26))
          S25[,,i] <- contrib
        }
      }
      S24[,,t] <- apply(S25, c(1,2), sum, na.rm = TRUE)
    }
    S27 <- apply(S24, c(1,2), sum, na.rm = TRUE)
    S28 <- S23 + S27
    S29 <- array(0, c(DimCovSusReInf, DimCovSusReInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      S30 <- array(0, c(DimCovSusReInf, DimCovSusReInf, NTotalpost))
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(ExpoTime[i] > t | ExpoTime[i] == 0){
          contrib <- -Pop[i] * CovSusReInf[GridIndic,] %*% t(CovSusReInf[GridIndic,]) *
            as.numeric(F2(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                 BetaCovInf, EstBetaCovSus, BetaCovSusReInf, alphaT) *
                         Av1(Rnd, GridIndic))
          S30[,,i] <- contrib
        }
      }
      S29[,,t] <- apply(S30, c(1,2), sum)
    }
    S31 <- apply(S29, c(1,2), sum, na.rm = TRUE)
    S32 <- array(0, c(DimCovSusReInf, DimCovSusReInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      S33 <- array(0, c(DimCovSusReInf, DimCovSusReInf, NTotalpost))
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          S34 <- sapply(1:NIterMC, function(L) Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                                                       i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L))
          contrib <- Pop[i] * CovSusReInf[GridIndic,] %*% t(CovSusReInf[GridIndic,]) *
            as.numeric(F2(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                 BetaCovInf, EstBetaCovSus, BetaCovSusReInf, alphaT) *
                         mean(S34))
          S33[,,i] <- contrib
        }
      }
      S32[,,t] <- apply(S33, c(1,2), sum, na.rm = TRUE)
    }
    S35 <- apply(S32, c(1,2), sum, na.rm = TRUE)
    S36 <- array(0, c(DimCovSusReInf, DimCovSusReInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      S37 <- array(0, c(DimCovSusReInf, DimCovSusReInf, NTotalpost))
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          S38 <- sapply(1:NIterMC, function(L) Av3(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                                                       i, GridIndic, t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L))
          contrib <- -Pop[i]^2 * CovSusReInf[GridIndic,] %*% t(CovSusReInf[GridIndic,]) *
            as.numeric((F2(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                  BetaCovInf, EstBetaCovSus, BetaCovSusReInf, alphaT))^2 * mean(S38))
          S37[,,i] <- contrib
        }
      }
      S36[,,t] <- apply(S37, c(1,2), sum, na.rm = TRUE)
    }
    S39 <- apply(S36, c(1,2), sum, na.rm = TRUE)
    S40 <- S31 + S35 + S39 + 0.1
    CheckEnoughCases(S40, "area-level reinfection susceptibility")
    EstBetaCovSusReInf <- BetaCovSusReInf - solve(S40) %*% S28
    contrib_time <- function(i, t) {
      GridIndic <- NLableGrid[i]
      if (ExpoTime[i] > t | ExpoTime[i] == 0) {
        return(-Pop[i] * as.numeric(
          F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
               BetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, alphaT) *
            Av1(Rnd, GridIndic)
        ))
      } else if (is_exposed(ExpoTime, IncPeriod, t, i)) {
        IA <- sapply(1:NIterMC, function(L)
          Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t,
                BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)
        )
        return(Pop[i] * as.numeric(
          F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
               BetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, alphaT) * mean(IA)
        ))
      } else return(0)
    }
    T1 <- sapply(1:MaxTimePand, function(t)
      sum(sapply(1:NTotalpost, contrib_time, t = t), na.rm = TRUE)
    )
    # T1 already merges both branches, so it is the complete score; adding it to
    # itself double-counted it.
    T2 <- sum(T1, na.rm = TRUE)
    T4 <- rep(0, MaxTimePand)
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          T4[t] <- T4[t] - Pop[i]^2 * as.numeric(
            (F1(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                  BetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, alphaT))^2 *
              mean(sapply(1:NIterMC, function(L) Av3(NLableGrid, Rnd, Dist, alphaS, delta,
                                                       lambda1, i, GridIndic, t,
                                                       BetaCovInf, BetaCovSus, BetaCovSusReInf,
                                                       alphaT, L)))
          )
        }
      }
    }
    T5 <- sum(T4, na.rm = TRUE)
    T6 <- T2 + T5
    EstAlphaT <- alphaT - NewtonStep(T2, T6)
    T7 <- array(0, c(DimCovInf, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(ExpoTime[i] > t | ExpoTime[i] == 0){
          T7[,,t] <- T7[,,t] - Pop[i] * F3(NLableGrid, Dist, EstAlphaS, delta,
                                             i, GridIndic, t, BetaCovInf,
                                             EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
            Av1(Rnd, GridIndic)
        }
      }
    }
    T8 <- apply(T7, c(1,2), sum, na.rm = TRUE)
    T9 <- array(0, c(DimCovInf, 1, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          T9[,,t] <- T9[,,t] + Pop[i] * F3(NLableGrid, Dist, EstAlphaS, delta,
                                             i, GridIndic, t, BetaCovInf,
                                             EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
            mean(sapply(1:NIterMC, function(L)
              Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                    i, GridIndic, t, BetaCovInf, BetaCovSus,
                    BetaCovSusReInf, alphaT, L)
            ))
        }
      }
    }
    T10 <- apply(T9, c(1,2), sum, na.rm = TRUE)
    T11 <- T8 + T10
    T12 <- array(0, c(DimCovInf, DimCovInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(ExpoTime[i] > t | ExpoTime[i] == 0){
          T12[,,t] <- T12[,,t] - Pop[i] * F4(NLableGrid, Dist, EstAlphaS, delta,
                                                   i, GridIndic, t, BetaCovInf,
                                                   EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
            as.numeric(Av1(Rnd, GridIndic))
        }
      }
    }
    T13 <- apply(T12, c(1,2), sum, na.rm = TRUE)
    T14 <- array(0, c(DimCovInf, DimCovInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          T14[,,t] <- T14[,,t] + Pop[i] * F4(NLableGrid, Dist, EstAlphaS, delta,
                                                   i, GridIndic, t, BetaCovInf,
                                                   EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
            mean(sapply(1:NIterMC, function(L)
              Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                    i, GridIndic, t, BetaCovInf, BetaCovSus,
                    BetaCovSusReInf, alphaT, L)
            ))
        }
      }
    }
    T15 <- apply(T14, c(1,2), sum, na.rm = TRUE)
    T16 <- T13 + T15
    T17 <- array(0, c(DimCovInf, DimCovInf, MaxTimePand))
    for(t in 1:MaxTimePand){
      for(i in 1:NTotalpost){
        GridIndic <- NLableGrid[i]
        if(is_exposed(ExpoTime, IncPeriod, t, i)){
          T17[,,t] <- T17[,,t] - Pop[i]^2 *
            F3(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t, BetaCovInf,
                 EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) %*%
            t(F3(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t, BetaCovInf,
                   EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT)) *
            mean(sapply(1:NIterMC, function(L)
              Av3(NLableGrid, Rnd, Dist, alphaS, delta, lambda1,
                    i, GridIndic, t, BetaCovInf, BetaCovSus,
                    BetaCovSusReInf, alphaT, L)
            ))
        }
      }
    }
    T18 <- apply(T17, c(1,2), sum, na.rm = TRUE)
    T19 <- T16 + T18
    CheckEnoughCases(T19, "individual-level infectivity")
    EstBetaCovInf <- BetaCovInf - solve(T19) %*% T11
    T20 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        GridIndic <- NLableGrid[i]
        if (ExpoTime[i] > t | ExpoTime[i] == 0) {
          Pop[i] * as.numeric(
            F5(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                     EstBetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
              Av1(Rnd, GridIndic)
          )
        } else 0
      }), na.rm = TRUE)
    }))
    T21 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        GridIndic <- NLableGrid[i]
        if (is_exposed(ExpoTime, IncPeriod, t, i)) {
          -Pop[i] * as.numeric(
            F5(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                     EstBetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
              mean(sapply(1:NIterMC, function(L)
                Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic,
                      t, BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)
              ))
          )
        } else 0
      }))
    }))
    T22 <- T20 + T21
    T23 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        GridIndic <- NLableGrid[i]
        if (ExpoTime[i] > t | ExpoTime[i] == 0) {
          -Pop[i] * as.numeric(
            F6(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                     EstBetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
              Av1(Rnd, GridIndic)
          )
        } else 0
      }))
    }))
    T24 <- rep(0, MaxTimePand)
    for(t in 1:MaxTimePand){
      T24[t] <- sum(sapply(1:NTotalpost, function(i){
        sum(sapply(1:NTotalGrid, function(GridIndic){
          if(NLableGrid[i] == GridIndic && ExpoTime[i] <= t && (ExpoTime[i] + IncPeriod[i]) > t && ExpoTime[i] != 0){
            Pop[i]*as.numeric(F6(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                       EstBetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT) *
                                mean(sapply(1:NIterMC, function(L)
                                  Av2(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t,
                                        BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L)))) -
              Pop[i]^2*as.numeric((F5(NLableGrid, Dist, EstAlphaS, delta, i, GridIndic, t,
                                            EstBetaCovInf, EstBetaCovSus, EstBetaCovSusReInf, EstAlphaT))^2 *
                                    mean(sapply(1:NIterMC, function(L)
                                      Av3(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, i, GridIndic, t,
                                            BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT, L))))
          } else 0
        }))
      }))
    }
    T25 <- sum(T24, na.rm = TRUE)
    T26 <- T23 + T25
    Estdelta <- delta - NewtonStep(T22, T26)
    # tau^2 * (lambda*D + (1-lambda)*I) is the PRECISION matrix: it is inverted
    # before being used as a covariance when the random effects are simulated, so
    # the density has to be evaluated against its inverse too.
    LGL1 <- function(par) {
      lambda1 <- par[1]
      tau1 <- par[2]
      if (lambda1 <= 0 || lambda1 >= 1 || tau1 <= 0) return(1e10)
      Sigma <- try(solve(tau1^2 * (lambda1 * D + (1 - lambda1) * diag(NTotalGrid))), silent = TRUE)
      if (inherits(Sigma, "try-error")) return(1e10)
      val <- -mean(sapply(1:NIterMC, function(L) dmvnorm(Rnd[L,], rep(0, NTotalGrid), sigma = Sigma, log = TRUE)))
      if (is.finite(val)) val else 1e10
    }
    EstU1fit <- optim(c(lambda1, tau1), fn = LGL1, method = "L-BFGS-B",
                      lower = c(1e-6, 1e-6), upper = c(1 - 1e-6, Inf))
    if (EstU1fit$convergence != 0) warning("The spatial variance components did not converge.", call. = FALSE)
    EstU1 <- EstU1fit$par
    EstGammau <- EstU1[1]
    HatSigmmaU <- EstU1[2]
    result=list(Rnd=Rnd,BetaCovInf=EstBetaCovInf,BetaCovSus=EstBetaCovSus,BetaCovSusReInf=EstBetaCovSusReInf,Uhat=EstU1,alphaS=EstAlphaS,alphaT=EstAlphaT,delta=Estdelta,tau1=HatSigmmaU,lambda1=EstGammau)
    result
  }
  O1=numeric()
  LGL1 <- function(NLableGrid, rndeft, Dist, alphaS, delta, lambda1, tau1,
                     BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) {
    O2 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        sum(sapply(1:NTotalGrid, function(GridIndic) {
          if (NLableGrid[i] == GridIndic && (ExpoTime[i] > t || ExpoTime[i] == 0)) {
            d1 <- mean(exp(rndeft[, GridIndic]))
            -Pop[i] * F1(NLableGrid, Dist, alphaS, delta, i, GridIndic, t,
                           BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT) * d1
          } else 0
        }))
      }))
    }))
    O3 <- sum(sapply(1:MaxTimePand, function(t) {
      sum(sapply(1:NTotalpost, function(i) {
        sum(sapply(1:NTotalGrid, function(GridIndic) {
          if (NLableGrid[i] == GridIndic &&
              ExpoTime[i] <= t && (ExpoTime[i] + IncPeriod[i]) > t && ExpoTime[i] != 0) {
            mean(sapply(1:NIterMC, function(L) {
              contrib <- sapply(1:NTotalpost, function(j) {
                if (NewLabelGrid[j, GridIndic] != 0 &&
                    InfTime[j] <= t && (InfTime[j] + InfPeriod[j]) >= t && InfTime[j] != 0) {
                  NInf[j] * exp(alphaT + CovInf[j, ] %*% BetaCovInf) * Dist[i, j]^(-delta)
                } else 0
              })
              dx <- sum(contrib[is.finite(contrib)])
              q <- exp(-Pop[i] * exp(alphaS + CovSus[GridIndic, ] %*% BetaCovSus +
                                       CovSusReInf[GridIndic, ] %*% BetaCovSusReInf +
                                       rndeft[L, GridIndic]) * dx)
              if (!is.finite(q) || q >= 1) 0 else log1p(-q)
            }))
          } else 0
        }))
      }))
    }))
    SigmaU <- solve(tau1^2 * (lambda1 * D + (1 - lambda1) * diag(NTotalGrid)))
    O4 <- mean(sapply(1:NIterMC, function(L) {
      dmvnorm(rndeft[L, ], rep(0, NTotalGrid), sigma = SigmaU, log = TRUE)
    }))
    O4 + O3 + O2
  }
  est0=estfun(NLableGrid,Dist,alphaS0,delta0,lambda0,tau0,BetaCovInf0,BetaCovSus0,BetaCovSusReInf0,alphaT0)
  alphaS=est0$alphaS
  delta=est0$delta
  lambda1=est0$lambda1
  tau1=est0$tau1
  BetaCovInf=est0$BetaCovInf
  BetaCovSus=est0$BetaCovSus
  BetaCovSusReInf=est0$BetaCovSusReInf
  alphaT=est0$alphaT
  Uhat=est0$Uhat
  rndeft=est0$Rnd
  AIC <- numeric()
  mes <- numeric()
  tolerance <- 0.1
  for (crtr in 1:NIterMCECM) {
    est <- estfun(NLableGrid, Dist, alphaS, delta, lambda1, tau1,
                  BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)
    list2env(est, envir = environment())
    O1[crtr] <- LGL1(NLableGrid, Rnd, Dist, alphaS, delta, lambda1, tau1,
                     BetaCovInf, BetaCovSus, BetaCovSusReInf, alphaT)
    # k = alphaS, alphaT, delta, tau, lambda plus the three beta vectors.
    NPar <- 5 + DimCovInf + DimCovSus + DimCovSusReInf
    AIC[crtr] <- -2 * O1[crtr] + 2 * NPar
    current_params <- c(alphaS, BetaCovInf, BetaCovSus, BetaCovSusReInf,
                        alphaT, delta, tau1, lambda1)
    mes <- if(crtr > 1) sqrt(sum((current_params - prev_params)^2)) else Inf
    prev_params <- current_params
    if(mes < tolerance) {
      message("MCECM Converged at iteration ", crtr, " with parameter change: ", mes)
      break
    }
    if(crtr %% 10 == 0) message("Iteration ", crtr, ": Parameter change = ", mes)
  }
  out1 <- list(
    alphaS = alphaS,
    BetaCovInf = BetaCovInf,
    BetaCovSus = BetaCovSus,
    BetaCovSusReInf = BetaCovSusReInf,
    alphaT = alphaT,
    delta = delta,
    tau1 = tau1,
    lambda1 = lambda1,
    AIC = AIC[crtr]
  )
  out1
}