R/population_simulation.R

In mortAAR: Analysis of Archaeological Mortality Data

#' Simulation of a population of adults with Gompertz distribution
#'
#' In many instances, it is useful to calculate with a population with
#' known parameters. To generate a population with realistic
#' characteristics is less obvious than it seems. We operate here
#' with the Gompertz distribution which provides a reasonable
#' approximation of human mortality for \bold{adult} mortality, that is
#' for the ages >= 15 years. The user has to specify
#' either the parameter b or the modal age M. The modal age M is
#' particular useful as it provides an intuitive understanding of
#' the resulting age distribution. In both instances, the second
#' parameter a is generated by the regression formula found by
#' \emph{Sasaki and Kondo 2016}. If neither is given, a population
#' with random parameters realistic for pre-modern times is generated.
#'
#' @param n number of individuals to be simulated.
#'
#' @param b numeric, optional. Gompertz parameter controlling the
#' level of mortality.
#'
#' @param M numeric, optional. Modal age M.
#'
#' @param start_age numeric. Start age, default: 15 years.
#'
#' @param max_age numeric. Maximal age, to avoid unlikely centenaries,
#' default: 100 years.
#'
#' @return
#' A list of two data.frames with the following items:
#'
#' First data.frame
#' \itemize{
#'   \item \bold{N}:  Number of individuals.
#'   \item \bold{b}:  Gompertz parameter controlling mortality.
#'   \item \bold{M}:  Modal age.
#'   \item \bold{a}:  Gompertz parameter controlling hazard of the
#'   youngest age group.
#'  }
#'
#'Second data.frame
#' \itemize{
#'   \item \bold{ind}:  ID of individuals.
#'   \item \bold{age}:  Simulated absolute age.
#'  }
#'
#' @references
#'
#' \insertRef{sasaki_kondo_2016}{mortAAR}
#'
#' @examples
#'
#' pop_sim <- pop.sim.gomp(n = 10000, M = 35)
#' pop_sim <- pop.sim.gomp(n = 10000, b = 0.03)
#' pop_sim <- pop.sim.gomp(n = 10000)

#' @rdname pop.sim.gomp
#' @export
pop.sim.gomp <- function(n, b = NULL, M = NULL, start_age = 15, max_age = 100) {
  if ( length(M) > 0) {
    M_1 <- 0
    M_2 <- 0
    while ( M < M_1 | M > M_2 ) {
      b <- stats::runif(n = 1, min = 0.02, max = 0.1)
      a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), 0.0823))
      M_ <- 1 / b * log (b/a) + start_age
      M_1 <- M_ - 0.001
      M_2 <- M_ + 0.001
    }
  } else if (length(b) > 0) {
    a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
    M <- 1 / b * log (b/a) + start_age
  } else {
    b <- stats::runif(n = 1, min = 0.02, max = 0.05)
    a <- exp(stats::rnorm(1, (-66.77 * (b - 0.0718) - 7.119), sqrt(0.0823) ) )
    M <- 1 / b * log (b/a) + start_age
  }
  count = 1
  lt_result <- matrix(NA, n, 2)
  while (count < (n + 1)) {
    age <- round(flexsurv::rgompertz(1, b, a) ) + start_age
    if (age < max_age) {
    lt_result[count,] <- c(ind = count, age = age)
    count <- count + 1
    }
  }
  lt_result <- data.frame(lt_result)
  colnames(lt_result) <- c("ind", "age")
  lt_values <- data.frame(n = n, b = b, a = a, M = M, start_age = start_age,
                          max_age = max_age)
  lt_list <- list(values = lt_values, result = lt_result)
  return(lt_list)
}