R/postNt.R

In MobilePhoneESSnetBigData/pestim: Population Estimations Using Mobile Phone Data

#' @title Posterior mean, median, and mode for the number of individuals at an arbitrary time.
#'
#' @description Compute the posterior mean, median, and mode for the number of individuals
#' generating posterior distribution according to the hierarchical model.
#'
#' @param nMNOmat transition matrix with the number of individuals displaced from cell to cell
#' detected by the Mobile Network Operator
#'
#' @param nReg non-negative integer vector with the number of individuals detected in each
#' cell according to the population register
#'
#' @param fu, fv named lists with the prior marginal distributions of the two-dimensional points
#' for the Monte Carlo integration
#'
#' @param flambda named list with the prior distribution of the lambda parameter
#'
#' @param distNames character vector with the names of the prior distributions for each cell
#'
#' @param variation list of lists whose components are parameters providing a measure of variation
#' of each prior distribution
#'
#' @param scale numeric vector with the scale to count the number of individuals. Default value is 1
#'
#' @param n number of points to generate in the posterior distribution for the computation. Default
#' value is 1e3
#'
#' @param relTol relative tolerance in the computation of the \code{\link{kummer}} function. Default
#' value is \code{1e-6}
#'
#' @param nSim number of two-dimensional points to generate to compute the integral. Default value
#' is \code{1e4}
#'
#' @param nStrata integer vector of length 2 with the number of strata in each dimension. Default
#' value is \code{c(1, 1e2)}
#'
#' @param verbose logical (default \code{FALSE}) to report progress of the computation
#'
#' @param nThreads number (default the number of all cores, including logical cores) to use for computation
#'
#' @param alpha the significance level for accuracy measures. Default value is 0.05
#'
#' @return \code{postNt} computes the posterior mean, median, and mode of the posterior distribution
#' for each cell at an arbitrary time \eqn{t}. The function returns a matrix with the estimates in
#' columns and the cells in rows.
#'
#' @details The prior distributions are specified as named lists where the first component of each
#' list must be the name of distribution ('unif', 'triang', 'degen', 'gamma') and the rest of
#' components must be named according to the name of the parameters of the random generator of the
#' corresponding distribution according to:
#'
#'   \itemize{
#'
#'     \item unif: \code{xMin}, \code{xMax} for the minimum, maximum of the sampled interval.
#'     \item degen: \code{x0} for the degenerate value of the random variable.
#'     \item triang: \code{xMin}, \code{xMax}, \code{xMode} for minimum, maximum and mode (see
#'     \code{\link{qtriang}}).
#'     \item gamma: \code{scale} and \code{shape} with the same meaning as in \code{\link{rgamma}}.
#'   }

#'
#' @seealso \code{\link{rNt}}, \code{\link{postN0}}, \code{\link{postNtcondN0}}
#'
#' @examples
#' ## First, the inputs:
#'
#' #The transition matrix of individuals detected by the MNO
#' nMNOmat <- rbind(c(10, 3, 4), c(5, 21, 3), c(3, 9, 18))
#'
#' # Population at the initial time of each cell according to the population register
#' nReg <- c(90, 130, 101)
#'
#' # List of priors for u
#' u0 <- rowSums(nMNOmat) / nReg
#' cv_u0 <- 0.15
#' fu <- lapply(u0, function(u){
#'  umin <- max(0, u - cv_u0 * u)
#'  umax <- min(1, u + cv_u0 * u)
#'  output <- list('unif', xMin = umin, xMax = umax)
#'  return(output)
#' })
#'
#' # List of priors for v
#' v0 <- nReg
#' cv_v0 <- 0.10
#' fv <- lapply(v0, function(u){
#'   umin <- max(0, u - cv_v0 * u)
#'   umax <- u + cv_v0 * u
#'   output <- list('unif', xMin = umin, xMax = umax)
#'   return(output)
#' })
#'
#' # List of priors for lambda
#' cv_lambda <- 0.6
#' alpha <- 1 / cv_lambda**2 - 1
#' flambda <- lapply(v0, function(v){list('gamma', shape = 1 + alpha, scale = v / alpha)})
#'
#' # Names and parameters of priors for the transition probabilities
#' distNames <- rep('unif', 3)
#' variation <- rep(list(list(cv = 0.20)), 3)
#'
#' # It takes a couple of minutes.
#' postNt(nMNOmat, nReg, fu, fv, flambda, distNames, variation)
#'
#' @include rNt.R
#' @include utils.R
#'
#' @import HDInterval
#'
#' @export
postNt <- function(nMNOmat, nReg, fu, fv, flambda, distNames, variation, scale = 1, n = 1e3,
                   relTol = 1e-6, nSim = 1e3, nStrata = c(1, 1e2), verbose = FALSE, nThreads = RcppParallel::defaultNumThreads(), alpha=0.05){

  Ntmat <- rNt(n, nMNOmat, nReg, fu, fv, flambda, distNames, variation,
               scale, relTol, nSim, nStrata, verbose, nThreads)
  postMean <- Ntmat[, round(mean(N)), by = c('cellID')]
  setnames(postMean, 'V1', 'value')
  postMean[, variable := 'postMean']


  postSD<-Ntmat[,round(sd(N),2), by = c('cellID')]
  setnames(postSD, 'V1', 'value')
  postSD[, variable := 'postSD']

  postCV<-Ntmat[, round(sd(N) / mean(N) * 100, 2), by = c('cellID')]
  setnames(postCV, 'V1', 'value')
  postCV[, variable := 'postCV']

  postMedian <- Ntmat[, round(median(N)), by = c('cellID')]
  setnames(postMedian, 'V1', 'value')
  postMedian[, variable := 'postMedian']

  postMedian_CILB<-Ntmat[, equalTailedInt(N, alpha)['lower'], by = c('cellID')]
  setnames(postMedian_CILB, 'V1', 'value')
  postMedian_CILB[, variable := 'postMedian_CILB']

  postMedian_CIUB<-Ntmat[, equalTailedInt(N, alpha)['upper'], by = c('cellID')]
  setnames(postMedian_CIUB, 'V1', 'value')
  postMedian_CIUB[, variable := 'postMedian_CIUB']

  postMedianQuantileCV<-Ntmat[, round( IQR(N) / median(N) * 100, 2), by = c('cellID')]
  setnames(postMedianQuantileCV, 'V1', 'value')
  postMedianQuantileCV[, variable := 'postMedianQuantileCV']

  postMode <- Ntmat[, Mode(N), by = c('cellID')]
  setnames(postMode, 'V1', 'value')
  postMode[, variable := 'postMode']

  postMode_CILB<-Ntmat[, hdi(N, 1-alpha)['lower'], by = c('cellID')]
  setnames(postMode_CILB, 'V1', 'value')
  postMode_CILB[, variable := 'postMode_CILB']

  postMode_CIUB<-Ntmat[, hdi(N, 1-alpha)['upper'], by = c('cellID')]
  setnames(postMode_CIUB, 'V1', 'value')
  postMode_CIUB[, variable := 'postMode_CIUB']

  postModeQuantileCV<-Ntmat[, round( IQR(N) / Mode(N) * 100, 2), by = c('cellID')]
  setnames(postModeQuantileCV, 'V1', 'value')
  postModeQuantileCV[, variable := 'postModeQuantileCV']


  DT <- rbindlist(list(postMean, postSD, postCV, postMedian, postMedian_CILB, postMedian_CIUB, postMedianQuantileCV, postMode,postMode_CILB, postMode_CIUB, postModeQuantileCV))
  output <- dcast(DT, cellID ~ variable, value.var = 'value')
  output[, cellID := NULL]
  setcolorder(output, c('postMean', 'postSD', 'postCV', 'postMedian', 'postMedian_CILB', 'postMedian_CIUB', 'postMedianQuantileCV', 'postMode', 'postMode_CILB', 'postMode_CIUB', 'postModeQuantileCV'))
  output <- as.matrix(output)[]
  return(output)

}