R/basta.R

Defines functions .CalcNonParamSurv .CalcLifeTable .CalcCVx .CalcGx .CalcHx .CalcEx .CalcDemoFunQuan .CalcKulbackLeibler .ExtractParalOut .RunMCMC .FillMCMCoutObj .CreateMCMCoutObj .SetParJumps .UpdateJumps .AssignMinAge.noMinAge .AssignMinAge.minAge .AssignMinAge .SampleAges.agecensus .SampleAges.agecmr .SampleAges .SamplePars .JitterPars .CalcHastRatio .CalcPriorAgeDist .SetPriorAgeDist .CalcPost .CreatePost .CalcLike.bastacensus .CalcLike.bastacmr .CalcLike .BuildLikeObj .DefineSurvNumeric .DefineSurvMatrix .DefineCumHazNumeric .DefineCumHazMatrix .DefineMortNumeric .DefineMortMatrix .ptnorm .dtnorm .rtnorm .CalcCovPars .CreateParObj .CreateFullParObj .SetDefaultTheta .CreateUserPar .MakeCovMat .FindCovType .CreateCovObj .BuildAliveMatrix .CreateAgeObj .CreateDataObj .CreateAlgObj coef.multibasta summary.multibasta print.multibasta multibasta CensusToCaptHist summary.basta print.basta plot.basta basta.default basta FixCMRdata print.bastaCheckCens print.bastaCheckCMR summary.bastaCheckCens summary.bastaCheckCMR DataCheck

Documented in basta basta.default CensusToCaptHist coef.multibasta DataCheck FixCMRdata multibasta plot.basta print.basta print.bastaCheckCens print.bastaCheckCMR print.multibasta summary.basta summary.bastaCheckCens summary.bastaCheckCMR summary.multibasta

# ============================== CODE METADATA =============================== #
# Functions to analyse mortality from capture-mark-recapture (CMR)
# or from census data 
# Name:        basta.R
# Date:        2022-10-05
# Version:     2.0
# Created by:  Fernando Colchero
# Modified by: Fernando Colchero

# General Comments: 
# ----------------- #
# - Combine packages BaSTA and BaSTA.ZIMS 
# - Extend the CMR part of the package to allow testing the effect of 
#   covariates in recapture probabilities (in progress).
# - Allow the CMR part of the package to have Min. and Max. Birth dates
# - Allow the CMR part to have times of censoring before the end of the 
#   study.

# List of names to change:
# ------------------------ #
# 1.- datObj to dataObj
# 2.- parObj$theta$priorMu to parObj$theta$priorMean
# 3.- parObj$gamma$priorMu to parObj$theta$priorMean
# 4.- parObj$eta$priorMu to parObj$lambda$priorMean

# ================================ CODE START ================================ #
# =========================== #
# === A) USER FUNCTIONS: ====
# =========================== #

# A.1) Data check function:
# ------------------------- #
# Data check function:
DataCheck <- function(object, dataType = "CMR", studyStart = NULL, 
                      studyEnd = NULL, silent = TRUE) {
  if (dataType == "CMR") {
    # Extract study duration:
    Ti <- studyStart
    Tf <- studyEnd
    st <- Ti:Tf
    nt <- length(st)
    
    # Individual IDs:
    idnames <- object[, 1]
    
    # Number of records:
    n <- nrow(object)
    
    # Birth and death:
    bd <- as.matrix(object[, 2:3])
    
    # Recapture matrix:
    Y <- as.matrix(object[, 1:nt + 3]); colnames(Y) <- st
    Y1 <- Y
    Y1[which(Y > 1)] <- 1
    Tm <- matrix(st, n, nt, byrow = TRUE)
    
    # Covariates:
    if(ncol(object) > nt + 3){
      Z <- as.matrix(object[, (nt + 4):ncol(object)])
    } else {
      Z <- matrix(1, n, 1)
    }
    
    # 1. Death before observations start
    type1 <- as.vector(which(bd[, 2] < Ti & bd[, 2] != 0 & !is.na(bd[, 2])))
    if (length(type1) != 0 & !silent) {
      cat("The following rows deaths occur before observations start:\n")
      print(type1)
    }
    
    # 2. No birth/death AND no obervations
    type2 <- as.vector(which(rowSums(bd) + rowSums(Y1) == 0))
    if (length(type2) != 0 & !silent) {
      cat("The following rows have no object (unknown birth, unknown death, and no observations):\n")
      print(type2)
    }
    
    # 3. Birth after death 
    type3 <- as.vector(which(bd[, 1] > bd[, 2] & bd[, 1] != 0 & 
                               !is.na(bd[, 1]) & bd[, 2] != 0 & 
                               !is.na(bd[, 2])))
    if (length(type3) != 0 & !silent) {
      cat("The following rows have birth dates that are later than their death dates:\n")
      print(type3)
    }
    
    # 4. Observations after death
    # Calculate first and last time observed: 
    st <- Ti:Tf
    ytemp <- t(t(Y1) * st)
    lastObs <- c(apply(ytemp, 1, max))
    tempDeath <- bd[, 2]
    tempDeath[which(tempDeath == 0)] <- Inf
    type4 <- as.vector(which(lastObs > tempDeath & tempDeath >= Ti))
    rm(tempDeath)
    
    if (length(type4) != 0 & !silent) {
      cat("The following rows have observations that occur after the year of death:\n")
      print(type4)
    }
    
    # 5. Observations before birth
    ytemp[ytemp == 0]<- Inf
    firstObs <- c(apply(ytemp, 1, min))
    type5 <- as.vector(which(firstObs < bd[, 1]))
    
    if (length(type5) != 0 & !silent) {
      cat("The following rows have observations that occur before the year of birth:\n")
      print(type5)
    }
    
    # 6. Year of birth should be a zero in recapture matrix Y
    idb <- which(bd[, 1] > 0 & !is.na(bd[, 1]) & bd[, 1] >= Ti & bd[, 1] <= Tf)
    bcol <- bd[idb, 1] - Ti
    bpos <- bcol * n + idb
    type6 <- as.vector(idb[which(Y[bpos] == 1)])
    
    if (length(type6) != 0 & !silent) {
      cat("The following rows have a one in the recapture matrix in the birth year:\n")
      print(type6)
    }
    
    n <- nrow(Y)   
    
    # All OK:
    if(length(c(type1, type2, type3, type4, type5, type6)) > 0) {
      stopExec <- TRUE
    } else {
      stopExec <- FALSE
    }
    if (!silent) {
      if (stopExec) {
        cat("Problems were detected with the data\nYou can use function FixCMRdata() to resolve issues.\n\n")
      } else {
        cat("No problems were detected with the data.\n\n")
      }
    }
    
    # Summary information:
    firstBirth <- NA
    lastBirth <- NA
    idb <- which(bd[, 1] > 0)
    if (length(idb) > 0) {
      firstBirth <- min(bd[idb, 1])
      lastBirth <- max(bd[idb, 1])
    }
    
    firstDeath <- NA
    lastDeath <- NA
    idd <- which(bd[, 2] > 0)
    if (length(idd) > 0) {
      firstDeath <- min(bd[idd, 2])
      lastDeath <- max(bd[idd, 2])
    }
    datSumm <- list(n = n, 
                    knownBirth = length(which(bd[, 1] > 0)),
                    knownDeath = length(which(bd[, 2] > 0)),
                    knownBD = length(which(bd[, 2] > 0 & bd[, 1] > 0)),
                    detects = sum(Y1), firstDetec = min(ytemp), 
                    lastDetec = max(ytemp[which(ytemp != Inf)]),
                    firstBirth = firstBirth,
                    lastBirth = lastBirth,
                    firstDeath = firstDeath,
                    lastDeath = lastDeath)
    
    if (ncol(object) > nt + 3) {
      newData <- data.frame(idnames, bd, Y, Z)
    } else {
      newData <- data.frame(idnames, bd, Y)
    }
    probDescr <- c(type1 = "Death < Entry", type2 = "No obs.", 
                   type3 = "Birth > Death", type4 = "Death < Depart",
                   type5 = "Birth > Entry", type6 = "Birth = 1 in Y")
    prbList <- list(newData = newData, type1 = type1, type2 = type2, 
                    type3 = type3, type4 = type4, type5 = type5, type6 = type6,
                    summary = datSumm, stopExec = stopExec, 
                    probDescr = probDescr, dataType = dataType, 
                    studyStart = studyStart, studyEnd = studyEnd)
    class(prbList) <- "bastaCheckCMR"
  } else if (dataType == "census") {
    prbList <- list()
    prbList$n <- nrow(object)
    prbn <- 0
    # Find relevant column:
    colsnames <- c("Birth.Date", "Min.Birth.Date", "Max.Birth.Date", 
                   "Entry.Date", "Depart.Date", "Depart.Type")
    idcols <- which(!colsnames %in% colnames(object))
    prbList$stopExec <- FALSE
    if (length(idcols) > 0) {
      prbn <- prbn + 1
      if (!silent) {
        cat("\nCritical column names missing:\n")
        for (i in 1:length(idcols)) {
          cat(sprintf("Column '%s' missing.\n", colsnames[idcols]))
        }
        cat("Column names required are:\n")
        cat(sprintf("%s.\n", paste(colsnames, collapse = ", ")))
      }
      prbList$missCols <- idcols
      prbList$stopExec <- TRUE
    } else {
      # Convert to dates:
      object$Birth.Date <- as.Date(object$Birth.Date, format = "%Y-%m-%d")
      object$Min.Birth.Date <- as.Date(object$Min.Birth.Date, format = "%Y-%m-%d")
      object$Max.Birth.Date <- as.Date(object$Max.Birth.Date, format = "%Y-%m-%d")
      object$Entry.Date <- as.Date(object$Entry.Date, format = "%Y-%m-%d")
      object$Depart.Date <- as.Date(object$Depart.Date, format = "%Y-%m-%d")
      
      # Find NA records in dates:
      prbList$nas <- list()
      nnas <- 0
      for (i in 1:length(colsnames)) {
        idna <- which(is.na(object[[colsnames[i]]]))
        if (length(idna) > 0) {
          prbList$nas[[colsnames[i]]] <- idna
          nnas <- nnas + 1
        } else {
          prbList$nas[[colsnames[i]]] <- "None"
        }
      }
      if (nnas > 0) {
        if (!silent) {
          cat(sprintf("\nThere are %s columns with NA records\n", nnas))
          cat("Use print(object) to review NA records.\n")
        }
        prbList$stopExec <- TRUE
      }
      
      # Find date ranges:
      prbList$DateRan <- apply(object[, colsnames[-6]], 2, range, na.rm = TRUE)
      
      # Find inconsistencies between dates:
      datepr <- c("Min Birth > Birth", "Birth > Max Birth", 
                  "Min Birth > Max Birth", "Birth > Entry", "Min Birth > Entry",
                  "Max Birth > Entry", "Entry > Depart")
      nchpr <- nchar(datepr)
      maxch <- max(nchpr)
      prbList$probDescr <- datepr
      prbList$MinBBirth <- which(object$Min.Birth.Date > object$Birth.Date)
      prbList$BirthMaxB <- which(object$Birth.Date > object$Max.Birth.Date)
      prbList$MinBMaxB <- which(object$Min.Birth.Date > object$Max.Birth.Date)
      prbList$BirthEntr <- which(object$Birth.Date > object$Entry.Date)
      prbList$MinBEntr <- which(object$Min.Birth.Date > object$Entry.Date)
      prbList$MaxBEntr <- which(object$Max.Birth.Date > object$Entry.Date)
      prbList$EntrDep <- which(object$Entry.Date > object$Depart.Date)
      idlens <- grep("MinBBirth", names(prbList)):grep("EntrDep", names(prbList))
      dateslen <- sapply(idlens, function(pr) length(prbList[[pr]]))
      names(dateslen) <- names(prbList)[idlens]
      if (any(dateslen > 0)) {
        if (!silent) {
          cat("\nThe following number of records\nhave inconsistencies between dates:\n")
          for (i in which(dateslen > 0)) {
            cat(sprintf("%s%s: %s records\n", 
                        datepr[i], paste(rep(" ", maxch + 2 - nchpr[i]),
                                         collapse = ""), dateslen[i]))
          }
          cat("\nUse print(object) to find inconsistent records.\n")
        }
        prbList$stopExec <- TRUE
      } else {
        if (!silent) {
          cat("\nNo inconsistencies between dates.\n")
        }
      }
      # Find types of departue types
      prbList$DepartType <- table(object$Depart.Type)
      departType <- as.character(object$Depart.Type)
      if (!all(departType %in% c("C", "D"))) {
        prbList$stopExec <- TRUE
      }
      
      # Uncensored individuals:
      prbList$idUnCens <- which(departType == "D")
      prbList$nUnCens <- length(prbList$idUnCens)
      
      # Censored individuals:
      prbList$idCens <- which(departType == "C")
      prbList$nCens <- length(prbList$idCens)
      
      # Birth times to be estimated:
      prbList$idNoBirth <- which(object$Birth.Date != object$Min.Birth.Date)
      prbList$nNoBirth <- length(prbList$idNoBirth)
    }
    
    # Output:
    class(prbList) <- "bastaCheckCens"
  } else {
    stop("Wrong 'dataType' specified. Options are 'CMR' and 'census'.")
  }
  return(prbList)
}

# Summary function for data check:
summary.bastaCheckCMR <- function(object, ...) {
  cat(sprintf("DATA CHECK: %s\n", Sys.time()))
  cat("================================\n")
  
  datSumm <- object$summary
  nchDsum <- sapply(1:length(datSumm), function(dd) {
    nch <- nchar(datSumm[[dd]])
    if (is.na(nch)) nch <- 2
    return(nch)
  })
  maxNchDsum <- max(nchDsum)
  summCats <- c("Number of individuals", "Number with known birth year",
                "Number with known death year", "Number with known birth\n  AND death years", "Number of recaptures", "Earliest detection year",
                "Latest detection year", "Earliest birth year",
                "Latest birth year", "Earliest death year", "Latest death year")
  nchCats <- nchar(summCats)
  nchCats[4] <- 15
  maxNchCats <- max(nchCats)
  tNch <- maxNchDsum + maxNchCats + 1
  
  cat("\nDATA SUMMARY:\n=============\n")
  for (ich in 1:length(datSumm)) {
    catLab <- sprintf("- %s:%s%s\n", summCats[ich], 
                      paste(rep(" ", tNch - nchCats[ich] - nchDsum[ich]), 
                            collapse = ""), datSumm[[ich]])
    cat(catLab)
  }
  
  cat("\nNUMBER OF PROBLEMATIC RECORDS:\n==============================\n")
  probDescr <- object$probDescr
  nchCats <- sapply(1:length(probDescr), function(dd) nchar(probDescr[dd]))
  nPrb <- sapply(1:length(probDescr), function(ipr) {
    ity <- sprintf("type%s", ipr)
    nipr <- length(object[[ity]])
    return(nipr)
  })
  nchPrb <- nchar(nPrb)
  tNch <- max(nchCats) + max(nchPrb) + 1
  
  if (object$stopExec) {
    for (ipr in 1:6) {
      catLab <- sprintf("- %s:%s%s\n", probDescr[ipr], 
                        paste(rep(" ", tNch - nchCats[ipr] - nchPrb[ipr]), 
                              collapse = ""), nPrb[[ipr]])
      cat(catLab)
    }
    cat("\nNote: you can use function FixCMRdata() to fix issues.\n\n")
  } else {
    cat("No problematic records detected.\n\n")
  }
}  

summary.bastaCheckCens <- function(object, ...) {
  cat(sprintf("DATA CHECK: %s\n", Sys.time()))
  cat("================================\n")
  
  if ("missCols" %in% names(object)) {
    cat("Data object is missing the following columns:\n")
    cat(sprintf("%s.\n", paste(object$missCols, collapse = ", ")))
    cat("The analysis cannot be performed without these columns.")
  } else {
    cat("\nDATA SUMMARY:\n=============\n")
    # Total number of records:
    cat(sprintf("Total number of records  : %s\n", object$n))
    # Number of censored records:
    cat(sprintf("Number censored (C)      : %s\n", object$nCens))
    # Number of censored records:
    cat(sprintf("Number uncensored (D)    : %s\n", object$nUnCens))
    # Number of records with unknown birth:
    cat(sprintf("Number with unknown birth: %s\n", object$nNoBirth))
    
    
    # Dates columns with NAs:
    cat("\nNAs IN DATES COLUMNS:\n")
    cat("---------------------\n")
    nalen <- sapply(1:length(object$nas), 
                    function(ii) {
                      if (object$nas[[ii]][1] != "None") {
                        nna <- length(object$nas[[ii]])
                      } else {
                        nna <- 0
                      }
                      return(nna)
                    })
    names(nalen) <- names(object$nas)
    idna <- which(nalen > 0)
    nachar <- nchar(names(object$nas))
    mnach <- max(nachar)
    if (any(nalen > 0)) {
      cat("NAs found in the following dates columns:\n")
      for (ii in idna) {
        cat(sprintf("%s:\n", names(object$nas)[ii]))
        print(object$nas[[ii]])
        cat("\n")
      }
    } else {
      cat("No NAs found in dates columns.\n")
    }
    
    
    # Dates ranges:
    cat("\nDATES RANGES:\n")
    cat("-------------\n")
    dtab <- object$DateRan
    idna <- which(is.na(dtab))
    if (length(idna) > 0) dtab[idna] <- "NANA-NA-NA"
    cat(sprintf("%s\n", paste(colnames(dtab), collapse = "\t")))
    for(ii in 1:2) {
      cat(sprintf("%s\n", paste(dtab[ii, ], collapse = "\t")))
    }
    
    
    # Inconsistencies in the dates columns:
    cat("\nINCONSISTENCIES BETWEEN DATES COLUMNS:\n")
    cat("--------------------------------------\n")
    
    datepr <- c("Min Birth > Birth", "Birth > Max Birth", 
                "Min Birth > Max Birth", "Birth > Entry", "Min Birth > Entry",
                "Max Birth > Entry", "Entry > Depart")
    idlens <- c("MinBBirth", "BirthMaxB", "MinBMaxB", "BirthEntr", "MinBEntr",
                "MaxBEntr", "EntrDep")
    nchpr <- nchar(datepr)
    maxch <- max(nchpr)
    dateslen <- sapply(idlens, function(pr) length(object[[pr]]))
    names(dateslen) <- names(object)[idlens]
    if (any(dateslen > 0)) {
      cat("Records with inconsistencies between dates:\n")
      for (i in which(dateslen > 0)) {
        cat(sprintf("\n%s:\n", datepr[i]))
        print(object[[idlens[i]]])
      }
    } else {
      cat("None.\n")
    }
    
    # Inconsistencies in Depart Type:
    cat("\nINCONSISTENCIES IN DEPARTURE TYPES:\n")
    cat("-----------------------------------\n")
    if (length(which(!names(object$DepartType) %in% c("C", "D")))) {
      cat(sprintf("Departure type codes found: %s\n", 
                  paste(names(object$DepartType), collapse = ", ")))
      cat("\nWarning: only departure types required\nare C (censored) and D (dead).")
    } else {
      cat("None\n")
    }
  }
  
}  

# Print function for data check:
print.bastaCheckCMR <- function(x, ...) {
  cat(sprintf("DATA CHECK: %s\n", Sys.time()))
  cat("================================\n")
  
  datSumm <- x$summary
  nchDsum <- sapply(1:length(datSumm), function(dd) {
    nch <- nchar(datSumm[[dd]])
    if (is.na(nch)) nch <- 2
    return(nch)
  })
  maxNchDsum <- max(nchDsum)
  summCats <- c("Number of individuals", "Number with known birth year",
                "Number with known death year", "Number with known birth\n  AND death years", "Number of recaptures", "Earliest detection year",
                "Latest detection year", "Earliest birth year",
                "Latest birth year", "Earliest death year", "Latest death year")
  nchCats <- nchar(summCats)
  nchCats[4] <- 15
  maxNchCats <- max(nchCats)
  tNch <- maxNchDsum + maxNchCats + 1
  
  cat("\nDATA SUMMARY:\n=============\n")
  for (ich in 1:length(datSumm)) {
    catLab <- sprintf("- %s:%s%s\n", summCats[ich], 
                      paste(rep(" ", tNch - nchCats[ich] - nchDsum[ich]), 
                            collapse = ""), datSumm[[ich]])
    cat(catLab)
  }
  
  cat("\nNUMBER OF PROBLEMATIC RECORDS:\n==============================\n")
  probDescr <- x$probDescr
  nchCats <- sapply(1:length(probDescr), function(dd) nchar(probDescr[dd]))
  nPrb <- sapply(1:length(probDescr), function(ipr) {
    ity <- sprintf("type%s", ipr)
    nipr <- length(x[[ity]])
    return(nipr)
  })
  nchPrb <- nchar(nPrb)
  tNch <- max(nchCats) + max(nchPrb) + 1
  
  if (x$stopExec) {
    for (ipr in 1:6) {
      catLab <- sprintf("- %s:%s%s\n", probDescr[ipr], 
                        paste(rep(" ", tNch - nchCats[ipr] - nchPrb[ipr]), 
                              collapse = ""), nPrb[[ipr]])
      cat(catLab)
    }
    cat("\nNote: you can use function FixCMRdata() to fix issues.\n\n")
  } else {
    cat("No problematic records detected.\n\n")
  }
}  

print.bastaCheckCens <- function(x, ...) {
  cat(sprintf("DATA CHECK: %s\n", Sys.time()))
  cat("================================\n")
  
  if ("missCols" %in% names(x)) {
    cat("Data object is missing the following columns:\n")
    cat(sprintf("%s.\n", paste(x$missCols, collapse = ", ")))
    cat("The analysis cannot be performed without these columns.")
  } else {
    cat("\nDATA SUMMARY:\n=============\n")
    # Total number of records:
    cat(sprintf("Total number of records  : %s\n", x$n))
    # Number of censored records:
    cat(sprintf("Number censored (C)      : %s\n", x$nCens))
    # Number of censored records:
    cat(sprintf("Number uncensored (D)    : %s\n", x$nUnCens))
    # Number of records with unknown birth:
    cat(sprintf("Number with unknown birth: %s\n", x$nNoBirth))
    
    
    # Dates columns with NAs:
    cat("\nNAs IN DATES COLUMNS:\n")
    cat("---------------------\n")
    nalen <- sapply(1:length(x$nas), 
                    function(ii) {
                      if (x$nas[[ii]][1] != "None") {
                        nna <- length(x$nas[[ii]])
                      } else {
                        nna <- 0
                      }
                      return(nna)
                    })
    names(nalen) <- names(x$nas)
    idna <- which(nalen > 0)
    nachar <- nchar(names(x$nas))
    mnach <- max(nachar)
    if (any(nalen > 0)) {
      cat("NAs found in the following dates columns:\n")
      for (ii in idna) {
        cat(sprintf("%s:\n", names(x$nas)[ii]))
        print(x$nas[[ii]])
        cat("\n")
      }
    } else {
      cat("No NAs found in dates columns.\n")
    }
    
    
    # Dates ranges:
    cat("\nDATES RANGES:\n")
    cat("-------------\n")
    dtab <- x$DateRan
    idna <- which(is.na(dtab))
    if (length(idna) > 0) dtab[idna] <- "NANA-NA-NA"
    cat(sprintf("%s\n", paste(colnames(dtab), collapse = "\t")))
    for(ii in 1:2) {
      cat(sprintf("%s\n", paste(dtab[ii, ], collapse = "\t")))
    }
    
    
    # Inconsistencies in the dates columns:
    cat("\nINCONSISTENCIES BETWEEN DATES COLUMNS:\n")
    cat("--------------------------------------\n")
    
    datepr <- c("Min Birth > Birth", "Birth > Max Birth", 
                "Min Birth > Max Birth", "Birth > Entry", "Min Birth > Entry",
                "Max Birth > Entry", "Entry > Depart")
    idlens <- c("MinBBirth", "BirthMaxB", "MinBMaxB", "BirthEntr", "MinBEntr",
                "MaxBEntr", "EntrDep")
    nchpr <- nchar(datepr)
    maxch <- max(nchpr)
    dateslen <- sapply(idlens, function(pr) length(x[[pr]]))
    names(dateslen) <- names(x)[idlens]
    if (any(dateslen > 0)) {
      cat("Records with inconsistencies between dates:\n")
      for (i in which(dateslen > 0)) {
        cat(sprintf("\n%s:\n", datepr[i]))
        print(x[[idlens[i]]])
      }
    } else {
      cat("None.\n")
    }
    
    # Inconsistencies in Depart Type:
    cat("\nINCONSISTENCIES IN DEPARTURE TYPES:\n")
    cat("-----------------------------------\n")
    if (length(which(!names(x$DepartType) %in% c("C", "D")))) {
      cat(sprintf("Departure type codes found: %s\n", 
                  paste(names(x$DepartType), collapse = ", ")))
      cat("\nWarning: only departure types required\nare C (censored) and D (dead).")
    } else {
      cat("None\n")
    }
  }
  
}  

# Function to fix CMR data issues:
FixCMRdata <- function(object, studyStart, studyEnd, autofix = rep(0, 6), 
                       silent = TRUE) {
  dcheck <- DataCheck(object = object, studyStart = studyStart, 
                      studyEnd = studyEnd, silent = TRUE)
  
  # Extract study duration:
  Ti <- studyStart
  Tf <- studyEnd
  st <- Ti:Tf
  nt <- length(st)
  
  # Individual IDs:
  idnames <- object[, 1]
  
  # Number of records:
  n <- nrow(object)
  
  # Birth and death:
  bd <- as.matrix(object[, 2:3])
  
  # Recapture matrix:
  Y <- as.matrix(object[, 1:nt + 3]); colnames(Y) <- st
  Y1 <- Y
  Y1[which(Y > 1)] <- 1
  Tm <- matrix(st, n, nt, byrow = TRUE)
  
  # Covariates:
  if(ncol(object) > nt + 3){
    Z <- as.matrix(object[, (nt + 4):ncol(object)])
  } else {
    Z <- matrix(1, n, 1)
  }
  
  # 1. Death before observations start
  type1 <- dcheck$type1
  if (length(type1) != 0) {
    cat("The following rows deaths occur before observations start:\n")
    print(type1)
    
    # Actions - remove those rows from bd, Y and Z
    if (autofix[1] == 1) {
      bd <- bd[-type1, ]
      Y <- Y[-type1, ]
      Y1 <- Y1[-type1, ]
      idnames <- idnames[-type1]
      Z <- Z[-type1, ]
      Tm <- Tm[-type1, ]
      n <- nrow(Y)
      cat("These records have been removed from the Dataframe\n")
    }
  }
  
  # 2. No birth/death AND no obervations
  type2 <- dcheck$type2
  if (length(type2) != 0) {
    cat("The following rows have no object (unknown birth, unknown death, and no observations):\n")
    print(type2)
    
    #Actions - remove those rows from bd, Y and Z
    if (autofix[2] == 1) {
      bd <- bd[-type2, ]
      Y <- Y[-type2, ]
      Y1 <- Y1[-type2, ]
      idnames <- idnames[-type2]
      Z <- Z[-type2, ]
      Tm <- Tm[-type2, ]
      n <- nrow(Y)
      cat("These records have been removed from the Dataframe\n")
    }
  }
  
  # 3. Birth after death 
  type3 <- dcheck$type3
  if (length(type3) != 0) {
    cat("The following rows have birth dates that are later than their death dates:\n")
    print(type3)
    
    # Actions - remove the death, birth, both records?
    if (autofix[3] == 1) {
      bd[type3,2] = 0; cat("The death records have been replaced with 0.\n\n")
    } else if (autofix[3] == 2) {
      bd[type3,1] = 0; cat("The birth records have been replaced with 0\n")
    } else if (autofix[3] == 3) {
      bd[type3,1:2] = 0; cat("The birth and death records have been replaced with 0\n")
    }
  }
  
  # 4. Observations after death
  # Calculate first and last time observed: 
  type4 <- dcheck$type4
  if (length(type4) != 0) {
    cat("The following rows have observations that occur after the year of death:\n")
    print(type4)
    
    # Actions - remove spurious post-death observations
    if (autofix[4] == 1) {
      Ymd <- ((Tm - bd[, 2]) * Y1)[type4, ]
      Ymd[Ymd > 0] <- 0
      Ymd[Ymd < 0] <- 1
      Y[type4,] <- Ymd
      Y1[type4,] <- Ymd
      cat("Observations that post-date year of death have been removed.\n\n")
    }
  }
  
  # 5. Observations before birth
  type5 <- dcheck$type5
  
  if (length(type5) != 0) {
    cat("The following rows have observations that occur before the year of birth:\n")
    print(type5)
    
    # Actions - remove spurious pre-birth observations
    if (autofix[5] == 1) {
      Ymd <- ((Tm - bd[, 1]) * Y1)[type5, ]
      Ymd[Ymd > 0] <- 1
      Ymd[Ymd < 0] <- 0
      Y[type5, ] <- Ymd
      Y1[type5, ] <- Ymd
      cat("Observations that pre-date year of birth have been removed.\n\n")
    }
  }
  
  # 6. Year of birth should be a zero in recapture matrix Y
  idb <- which(bd[, 1] > 0 & !is.na(bd[, 1]) & bd[, 1] >= Ti & bd[, 1] <= Tf)
  bcol <- bd[idb, 1] - Ti
  bpos <- bcol * n + idb
  type6 <- dcheck$type6
  if (length(type6) != 0) {
    cat("The following rows have a one in the recapture matrix in the birth year:\n")
    print(type6)
    
    # Actions - put a zero.
    if (autofix[6] == 1) {
      Y[bpos] <- 0
      Y1[bpos] <- 0
      cat("These cells have been changed to 0.\n\n")
    }
  }
  
  n <- nrow(Y)   
  
  # All OK:
  
  if (ncol(object) > nt + 3) {
    newData <- data.frame(idnames, bd, Y, Z)
  } else {
    newData <- data.frame(idnames, bd, Y)
  } 
  dcheck$newData <- newData
  dcheckNew <- DataCheck(object = newData, studyStart = studyStart, 
                         studyEnd = studyEnd)
  
  return(dcheckNew)
}

# A.2) main basta function:
# ------------------------- #
basta <- function(object, ... ) UseMethod("basta")

basta.default <- function(object, dataType = "CMR", 
                          model = "GO", shape = "simple", 
                          studyStart = NULL, studyEnd = NULL, minAge = 0,
                          covarsStruct = "fused", formulaMort = NULL, 
                          formulaRecap = NULL, recaptTrans = studyStart, 
                          niter = 22000, burnin = 2001, thinning = 40, 
                          nsim = 1, parallel = FALSE, ncpus = 2, 
                          updateJumps = TRUE, negSenescence = FALSE, ...) {
  # Create BaSTA environment:
  # bastaenv <- new.env()
  
  # Additional arguments:
  argList <- list(...)
  
  # Check dataset and stop if problems found:
  # dcheck <- DataCheck(object, silent = TRUE)
  dcheck <- list(stopExec = FALSE)
  if (dcheck$stopExec) {
    stop("\nProblems detected with the data.\n", 
         "Use DataCheck() to find problems.", call. = FALSE)
  }
  
  # Variables to pass to all cpus:
  intVars <- c("algObj", "dataObj", "ageObj", "covObj", "defTheta", 
               "fullParObj", "parObj", "parCovObj", 
               "updateJumps", "niter", "nsim", "burnin", "thinning", "minAge", 
               ".CalcMort", '.CalcSurv', ".CalcMort.numeric", 
               ".CalcMort.matrix", ".CalcSurv.numeric", ".CalcSurv.matrix", 
               ".CalcCumHaz", ".CalcCumHaz.matrix", ".CalcCumHaz.numeric")
  
  # ------------------- #
  # General data setup:
  # ------------------- #
  # Algorithm information:
  algObj <- .CreateAlgObj(dataType, model, shape, studyStart, studyEnd, 
                          minAge, covarsStruct, formulaMort, formulaRecap, 
                          recaptTrans, niter, burnin, thinning, updateJumps,
                          nsim, negSenescence)
  
  # Dataset object:
  dataObj <- .CreateDataObj(object, algObj)
  
  # Covariate object:
  covObj <- .CreateCovObj(object, dataObj, algObj)
  
  # Create covariate names object:
  covsNames <- list(cat = NA, con = NA, class = NA)
  if (inherits(covObj, c("cateCov", "bothCov"))) {
    covsNames$cat <- covObj$cat
  }
  if (inherits(covObj, c("contCov", "bothCov"))) {
    covsNames$con <- covObj$cont
  }
  covsNames$class <- class(covObj)[1]
  
  # Define theta parameters:
  defTheta <- .SetDefaultTheta(algObj)
  
  # Incorporate user parameters:
  userPars <- .CreateUserPar(covObj, argList)
  
  # Construct full parameter object:
  fullParObj <- .CreateFullParObj(covObj, defTheta, algObj, userPars, 
                                  dataObj)
  
  # Initial parameter object:
  parObj <- .CreateParObj(fullParObj)
  
  # Covariate-parameter matrices:
  parCovObj <- .CalcCovPars(parObj = parObj, parCov = NULL, 
                            covObj = covObj, dataObj = dataObj, 
                            type = "both")
  
  # Age object:
  ageObj <- .CreateAgeObj(dataObj, algObj)
  
  # indices of kept values:
  outidObj <- list(all = seq(1, algObj$niter, algObj$thinning))
  outidObj$nAll <- length(outidObj$all)
  outidObj$idKeep <- which(outidObj$all >= algObj$burnin)
  outidObj$keep <- outidObj$all[outidObj$idKeep]
  outidObj$nKeep <- length(outidObj$keep)
  
  # ---------------------- #
  # Demographic functions:
  # ---------------------- #
  # a) Mortality function:
  .CalcMort <- function(theta, ...) UseMethod(".CalcMort")
  .CalcMort.matrix <- .DefineMortMatrix(algObj)
  .CalcMort.numeric <- .DefineMortNumeric(algObj)
  
  # b) Cummulative hazard:
  .CalcCumHaz <- function(theta, ...) UseMethod(".CalcCumHaz")
  .CalcCumHaz.matrix <- .DefineCumHazMatrix(algObj)
  .CalcCumHaz.numeric <- .DefineCumHazNumeric(algObj)
  
  # c) Survival:
  .CalcSurv <- function(theta, ...) UseMethod(".CalcSurv")
  .CalcSurv.matrix <- .DefineSurvMatrix(algObj)
  .CalcSurv.numeric <- .DefineSurvNumeric(algObj)
  
  # --------------------------------------- #
  # Assign variables to global environment:
  # --------------------------------------- #
  # for (name in intVars) assign(name, get(name), envir = globalenv())
  
  # ------------------------------ #
  # Initial likelihood and priors:
  # ------------------------------ #
  # Likelihood:
  likeObj <- .CalcLike(dataObj, parCovObj, parObj, fullParObj, ageObj, 
                       likeObj = NULL, algObj, ind = 1:dataObj$n, 
                       .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                       .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                       .CalcCumHaz, .CalcCumHaz.numeric, .CalcCumHaz.matrix)
  
  # Priors age object:
  priorAgeObj <- .SetPriorAgeDist(fullParObj, dataObj, covObj, .CalcSurv,
                                  .CalcSurv.numeric, .CalcSurv.matrix)
  
  # include in intVars:
  intVars <- c(intVars, "priorAgeObj")
  # assign("priorAgeObj", get("priorAgeObj"), envir = globalenv())
  
  # ----------------------- #
  # Run MCMC to find jumps:
  # ----------------------- #
  cat("\nRunning sequence to find jump SDs... ")
  Start <- Sys.time()
  jumpRun <- .RunMCMC(1, UpdJumps = TRUE, parJumps = NA, ageObj, dataObj,
                      parObj, fullParObj, covObj, 
                      priorAgeObj, algObj, outidObj, 
                      .CalcMort, .CalcMort.numeric, 
                      .CalcMort.matrix, .CalcSurv, 
                      .CalcSurv.numeric, .CalcSurv.matrix,
                      .CalcCumHaz, .CalcCumHaz.numeric, 
                      .CalcCumHaz.matrix, .JitterPars)
  End <- Sys.time()
  cat("Done\n")
  compTime <- round(as.numeric(End-Start, units = units(End - Start)), 2)
  cat(sprintf("Total jump SDs computing time: %.2f %s.\n\n", compTime, 
              units(End - Start)))
  
  
  # -------------- #
  # Run main MCMC:
  # -------------- #
  Start <- Sys.time()
  if (nsim > 1) {
    cat("Multiple simulations started...\n\n") 
    if (parallel) {
      sfInit(parallel = TRUE, cpus = ncpus)
      sfExport(list = c(intVars, ".Random.seed"))
      # if (Sys.info()[4] == "ADM-130757-Mac") {
      #   sfSource(file = "/Users/colchero/FERNANDO/PROJECTS/4.PACKAGES/BaSTA2.0/tests/sourceBaSTA.R")
      # } else {
      #   sfSource(file = "/Users/fernando/FERNANDO/PROJECTS/4.PACKAGES/BaSTA2.0/tests/sourceBaSTA.R")
      # }
      sfLibrary("BaSTA", character.only = TRUE,
      warn.conflicts = FALSE)
      bastaOut <- sfClusterApplyLB(1:nsim, .RunMCMC, UpdJumps = FALSE, 
                                   parJumps = jumpRun$jumps, ageObj, dataObj,
                                   parObj, fullParObj, covObj, 
                                   priorAgeObj, algObj, outidObj, 
                                   .CalcMort, .CalcMort.numeric, 
                                   .CalcMort.matrix, .CalcSurv, 
                                   .CalcSurv.numeric, .CalcSurv.matrix,
                                   .CalcCumHaz, .CalcCumHaz.numeric, 
                                   .CalcCumHaz.matrix, .JitterPars)
      sfRemoveAll(hidden = TRUE)
      sfStop()
    } else {
      bastaOut <- lapply(1:nsim, .RunMCMC, UpdJumps = FALSE, 
                         parJumps = jumpRun$jumps, ageObj, dataObj,
                         parObj, fullParObj, covObj, priorAgeObj, 
                         algObj, outidObj, .CalcMort, .CalcMort.numeric,
                         .CalcMort.matrix, .CalcSurv, .CalcSurv.numeric,
                         .CalcSurv.matrix, .CalcCumHaz, .CalcCumHaz.numeric, 
                         .CalcCumHaz.matrix, .JitterPars)
    }
  } else {
    cat("Simulation started...\n\n")
    bastaOut <- lapply(1:nsim, .RunMCMC, UpdJumps = FALSE, 
                       parJumps = jumpRun$jumps, ageObj, dataObj, parObj, 
                       fullParObj, covObj, priorAgeObj, algObj, outidObj,
                       .CalcMort, .CalcMort.numeric,
                       .CalcMort.matrix, .CalcSurv, .CalcSurv.numeric,
                       .CalcSurv.matrix, .CalcCumHaz, .CalcCumHaz.numeric, 
                       .CalcCumHaz.matrix, .JitterPars)
  }
  End <- Sys.time()
  cat("Simulations finished.\n")
  compTime <- round(as.numeric(End-Start, units = units(End - Start)), 2)
  cat(sprintf("Total MCMC computing time: %.2f %s.\n\n", compTime, 
              units(End - Start)))
  names(bastaOut) <- paste("sim.", 1:nsim, sep = "")
  
  # Calculate summaries:
  bastaSumars <- .ExtractParalOut(bastaOut, outidObj, fullParObj, covObj, 
                                  covsNames, nsim, dataObj, algObj, defTheta, 
                                  .CalcMort, .CalcMort.numeric, 
                                  .CalcMort.matrix, .CalcSurv, 
                                  .CalcSurv.matrix, .CalcSurv.numeric)
  # Create final BaSTA object:
  bastaFinal <- bastaSumars
  
  # Include individual runs in output:
  bastaFinal$runs <- bastaOut
  
  # Summary information for parameters:
  bastaFinal$fullpar <- fullParObj
  bastaFinal$simthe <- defTheta
  bastaFinal$jumps <- jumpRun$jumps
  
  # Covariate information:
  bastaFinal$covs <- covsNames
  
  # MCMC settings:
  bastaFinal$settings <- c(niter, burnin, thinning, nsim)
  names(bastaFinal$settings) <- c("niter", "burnin", "thinning", "nsim")
  
  # Model specifications:
  bastaFinal$modelSpecs <- 
    c(model, shape, minAge, covarsStruct,
      paste(names(covObj$cat), collapse = ", "), 
      paste(names(covObj$cont), collapse = ", "), dataType)
  names(bastaFinal$modelSpecs) <- c("model", "shape", "min. age", 
                                    "Covar. structure", "Categorical", 
                                    "Continuous", "DataType")
  # Life table:
  bastaFinal$lifeTable <- .CalcLifeTable(bastaFinal, covObj, algObj,
                                         dataObj)
  
  # Remove variables from global environment:
  rm(list = intVars)
  
  # Assign class:
  class(bastaFinal) <- "basta"
  return(bastaFinal)
}

# A.3) plotting and printing BaSTA outputs:
# ----------------------------------------- #
# Plotting:
plot.basta <- function(x, type = "traces", trace.name = "theta",
                       noCIs = FALSE, minSurv = NULL, ...) {
  # Additional arguments:
  args <- list(...)
  namesArgs <- names(args)
  
  # -------------------------- #
  # ---- Legacy arguments ---- #
  # find whether plot.type was used:
  if ("plot.type" %in% namesArgs) {
    type <- args$plot.type
  }
  
  # Find whether densities were specified:
  if ("densities" %in% namesArgs) {
    type <- "densities"
  }
  # ---- End of legacy arguments ---- #
  # --------------------------------- #
  
  # User par settings:
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  
  # Wrong type:
  if (!type %in% c("traces", "densities", "demorates", "gof", "fancy")) {
    stop("Wrong value for argument 'type'.\nAlternatives are 'traces', 'demorates', 'densities', 'gof', or 'fancy'.")
  }
  
  nv <- ifelse(type == "traces", x$settings['nsim'], length(x$surv))
  
  if ("col" %in% namesArgs) {
    Palette <- args$col
    if (length(Palette) < nv) {
      warning("Insufficient number of colors. Not all traces will be displayed.",
              call. = FALSE)
    }
  } else {
    if (nv <= 9) {
      Palette <- c('#E41A1C', '#377EB8', '#4DAF4A', '#984EA3', '#FF7F00', 
                   '#FFFF33', '#A65628', '#F781BF', '#999999')
    } else {
      Palette <- rainbow(nv)
    }
  }
  if ("lwd" %in% namesArgs) {
    lwd <- args$lwd
  } else {
    lwd <- 1
  }
  if ("lty" %in% namesArgs) {
    lty <- args$lty
  } else {
    lty <- 1
  }
  
  if (type %in% c("demorates", "fancy")) {
    catNames <- names(x$mort)
    lenCat <- length(catNames)
    
    if ("names.legend" %in% namesArgs) {
      if (length(args$names.legend) != nv) {
        stop(sprintf("Wrong length in 'names.legend'. Correct length is %s elements.",
                     nv), call. = FALSE)
      } else {
        legNames <- args$names.legend
      }
    } else {
      legNames <- catNames
    }
  }
  
  # ========== #
  # PARAMETERS:
  # ========== #
  if (type %in% c("traces", "densities")) {
    nsim <- x$settings["nsim"]
    # ----------- #
    # densities:
    # ---------- #
    if (type == "densities") {
      if (trace.name == "theta") {
        if (x$covs$class %in% c("fused", "inMort")) {
          if (!is.na(x$covs$cat[1])) {
            ncat <- length(x$covs$cat)
          } else {
            ncat <- 1
          }
        } else {
          ncat <- 1
        }
        pcol <- 1
        idpars <- sapply(x$fullpar$theta$names, function(pn) {
          grep(pn, colnames(x$params))
        })
        npar <- x$simthe$length
        prow <- npar
      } else if (trace.name == "gamma") {
        idpars <- grep("gamma", colnames(x$params))
        if (inherits(x$fullpar, "theta")) {
          stop("'gamma' parameters not calculated (not propHaz).", 
               call. = FALSE)
        } else {
          npar <- x$fullpar$gamma$len
          pcol <- ceiling(npar / 2)
          prow <- ceiling(npar / pcol)
        }
      } else if (trace.name == "lambda") {
        if (inherits(x$fullpar, "noLambda")) {
          stop("'lambda' parameters not calculated (no minAge).", 
               call. = FALSE)
        } else {
          pcol <- 1
          prow <- 1
          idpars <- grep("lambda", colnames(x$params))
          npar <- 1
        }
      } else if (trace.name == "pi") {
        if (x$modelSpecs["DataType"] == "CMR") {
          pcol <- ceiling(npar / 2)
          prow <- ceiling(npar / pcol)
          npar <- x$fullpar$pi$len
          idpars <- grep("pi", colnames(x$params))
        } else {
          stop("Argument 'trace.name' cannot be 'pi'.\n", 
               "No recapture probabilities estimated, dataType is not CMR.",
               call. = FALSE)
        }
      } else {
        stop("Wrong 'trace.name' specified.\n", 
             "Names are 'theta', 'gamma', 'lambda', or 'pi'", call. = FALSE)
      }
      par (mfrow = c(prow, pcol), mar = c(4, 4, 1, 1))
      for (pp in 1:npar) {
        if (trace.name == "theta") {
          main <- x$simthe$name[pp]
        } else {
          main <- x$fullpar[[trace.name]]$names[pp]
        }
        idpp <- grep(main, colnames(x$params))
        ddlist <- list()
        xlim <- c(NA, NA)
        ylim <- c(0, NA)
        for (ppi in 1:length(idpp)) {
          ddlist[[ppi]] <- density(x$params[, idpp[ppi]])
          xlim[1] <- min(c(xlim[1], ddlist[[ppi]]$x), na.rm = TRUE)
          xlim[2] <- max(c(xlim[2], ddlist[[ppi]]$x), na.rm = TRUE)
          ylim[2] <- max(c(ylim[2], ddlist[[ppi]]$y), na.rm = TRUE)
        }
        plot(xlim, ylim, col = NA, xlab = "", ylab = "", 
             main = main)
        for (ppi in 1:length(idpp)) {
          xx <- ddlist[[ppi]]$x
          yy <- ddlist[[ppi]]$y
          CIs <- unlist(x$coefficients[idpp[ppi], c("Lower95%CI", "Upper95%CI")])
          idcis <- which(xx >= CIs[1] & xx <= CIs[2])
          lines(xx, yy, col = Palette[ppi])
          polygon(c(xx[idcis], rev(xx[idcis])), 
                  c(yy[idcis], rep(0, length(yy[idcis]))), 
                  col = adjustcolor(Palette[ppi], alpha.f = 0.25), border = NA)
          if (grepl("gamma", colnames(x$params)[idpp[ppi]])) {
            lines(c(0, 0), ylim, col = 'grey60', lty = 1)
          }
        }
      }
      if (x$covs$class %in% c("fused", "inMort") & trace.name == "theta") {
        legend("topright", legend = names(x$covs$cat), 
               col = Palette[1:length(idpp)],
               pch = 22, pt.bg = adjustcolor(Palette[1:length(idpp)], 
                                             alpha.f = 0.25), 
               pt.cex = 2, bty = "n")
      }
      
      # ------- #
      # traces:
      # ------- #
    } else {
      if (trace.name == "theta") {
        if (x$covs$class %in% c("fused", "inMort")) {
          if (!is.na(x$covs$cat[1])) {
            ncat <- length(x$covs$cat)
          } else {
            ncat <- 1
          }
        } else {
          ncat <- 1
        }
        pcol <- ncat
        prow <- ceiling(x$fullpar$theta$len / pcol)
        npar <- x$simthe$length
        allnpar <- x$fullpar$theta$len
      } else if (trace.name == "gamma") {
        npar <- ncol(x$runs$sim.1$gamma)
        if (is.null(npar)) {
          stop("'gamma' parameters not calculated (not propHaz)", call. = FALSE)
        }
        pcol <- ceiling(npar / 2)
        prow <- ceiling(npar / pcol)
        allnpar <- x$fullpar$gamma$len
      } else if (trace.name == "lambda") {
        npar <- ifelse("lambda" %in% colnames(x$params), 1, NA)
        if (is.na(npar)) {
          stop("'lambda' parameters not calculated (no minAge)", call. = FALSE)
        }
        pcol <- 1
        prow <- 1
        allnpar <- 1
      } else if (trace.name == "pi") {
        if (x$modelSpecs["DataType"] == "CMR") {
          npar <- ncol(x$runs$sim.1$pi)
          pcol <- ceiling(npar / 2)
          prow <- ceiling(npar / pcol)
          allnpar <- x$fullpar$pi$len
        } else {
          stop("Argument 'trace.name' cannot be 'pi'.\n", 
               "No recapture probabilities estimated, dataType is not CMR.",
               call. = FALSE)
        }
      } else {
        stop("Wrong 'trace.name' specified.\n", 
             "Names are 'theta', 'gamma', 'lambda', or 'pi'", call. = FALSE)
      }
      keep <- seq(1, x$setting["niter"], x$settings["thinning"])
      par (mfrow = c(prow, pcol), mar = c(4, 4, 1, 1))
      for (pp in 1:allnpar) {
        ylim <- 
          range(sapply(1:nsim, function(ii) 
            range(x$runs[[ii]][[trace.name]][, pp]))) 
        xlim <- c(0, x$setting["niter"])
        main <- x$fullpar[[trace.name]]$names[pp]
        plot(xlim, ylim, col = NA, xlab = "", ylab = "", 
             main = main)
        for (tt in 1:nsim) {
          lines(keep, x$runs[[tt]][[trace.name]][, pp], 
                col = Palette[tt])
        }
      }
    }
    # =========== #
    # DEMO RATES:
    # =========== #
  } else if (type == "demorates") {
    par(mfrow = c(2, 1), mar = c(4, 4, 1, 1)) 
    demvname <- c("Mortality", "Survival")
    names(demvname) <- c("mort", "surv")
    for (demv in c("mort", "surv")) {
      ylim <- c(0, 0)
      if ("xlim" %in% namesArgs) {
        xlim <- args$xlim
      } else {
        xlim <- c(0, 0)
      }
      minAge <- as.numeric(x$modelSpecs["min. age"])
      for (icat in 1:lenCat) {
        if (is.null(minSurv)) {
          cuts <- x$cuts[[icat]]
        } else {
          cuts <- which(x$surv[[icat]][1, ] >= minSurv)
        }
        ylim <- range(c(ylim, x[[demv]][[icat]][, cuts]), 
                      na.rm = TRUE)
        if (! "xlim" %in% namesArgs) {
          xlim <- range(c(xlim, x$x[cuts] + minAge), na.rm = TRUE)
        }
      }
      plot(xlim, ylim, col = NA, xlab = "", ylab = demvname[demv])
      if (minAge > 0) lines(rep(minAge, 2), ylim, lty = 2)
      nn <- 0
      for (icat in 1:lenCat) {
        nn <- nn + 1
        if (is.null(minSurv)) {
          cuts <- x$cuts[[icat]]
        } else {
          cuts <- which(x$surv[[icat]][1, ] >= minSurv)
        }
        yy <- x[[demv]][[icat]][, cuts]
        xx <- x$x[cuts]
        if (!noCIs) {
          polygon(c(xx, rev(xx)) + minAge, c(yy[2, ], rev(yy[3, ])), 
                  col = adjustcolor(Palette[nn], alpha.f = 0.25),
                  border = NA)
        }
        lwdd <- ifelse(length(lwd) > 1, lwd[nn], lwd)
        ltyy <- ifelse(length(lty) > 1, lty[nn], lty)
        lines(xx + minAge, yy[1, ], lwd = lwdd, col = Palette[nn], 
              lty = ltyy)
      }
      if (nn > 1 & demv == 'surv') {
        legend(x = 'topright', legend = legNames, 
               pch = 22, pt.cex = 2, cex = 1.25, 
               pt.bg = adjustcolor(Palette, alpha.f = 0.25), col = Palette,
               bty = 'n')
        
      }
    }
    # ================ #
    # GOODNESS OF FIT: 
    # ================ #
  } else if (type == "gof") {
    ncat <- length(x$surv)
    catname <- names(x$surv)
    pcol <- ceiling(ncat / 2)
    prow <- ceiling(ncat / pcol)
    pcol <- 2
    prow <- ncat
    
    par(mfrow = c(prow, pcol), mar = c(4, 4, 1, 1)) 
    for (icat in 1:ncat) {
      ylim <- c(0, 1)
      xlim <- c(0, 0)
      minAge <- as.numeric(x$modelSpecs["min. age"])
      if (is.null(minSurv)) {
        cuts <- x$cuts[[icat]]
      } else {
        cuts <- which(x$surv[[icat]][1, ] >= minSurv)
      }
      if (! "xlim" %in% namesArgs) {
        xlim <- range(c(xlim, x$x[cuts] + minAge), na.rm = TRUE)
      } else {
        xlim <- args$xlim
      }
      if (is.data.frame(x$lifeTable[[icat]])) {
        lifeTab <- x$lifeTable[[icat]]
        PLE <- FALSE
      } else {
        lifeTab <- x$lifeTable[[icat]]$Mean
        ple <- x$lifeTable[[icat]]$ple
        PLE <- TRUE
      }
      minAge <- as.numeric(x$modelSpecs["min. age"])
      idAges <- which(lifeTab$Ages >= minAge)
      ltMinAge <- lifeTab$Ages[idAges[1]]
      
      # Survival:
      if (catname[icat] == "nocov") {
        main <- ""
      } else {
        main <- catname[icat]
      }
      plot(xlim, ylim, col = NA, xlab = "Age", ylab = "Survival", 
           main = main)
      if (minAge > 0) lines(rep(minAge, 2), ylim, lty = 2, col = 'orange')
      nn <- 0
      yy <- x$surv[[icat]][, cuts]
      xx <- x$x[cuts]
      if (ltMinAge > minAge) {
        idxx <- which(xx + minAge >= ltMinAge)
        yy <- yy[, idxx] / yy[, idxx[1]]
        xx <- xx[idxx]
      }
      if (!noCIs) {
        polygon(c(xx, rev(xx)) + minAge, c(yy[2, ], rev(yy[3, ])), 
                col = adjustcolor(Palette[1], alpha.f = 0.25),
                border = NA)
      }
      lwdd <- ifelse(length(lwd) > 1, lwd[1], lwd)
      ltyy <- ifelse(length(lty) > 1, lty[1], lty)
      # ======= BUG 2023-10-19 ========= #
      lines(lifeTab$Ages[idAges], lifeTab$lx[idAges] / lifeTab$lx[idAges[1]],
            type = "s")
      # ================================ #
      
      if (PLE) {
        if (minAge > 0) {
          idpl <- which(ple$Ages >= ltMinAge)
          if (ple$Ages[idpl[1]] > ltMinAge) {
            pleAges <- c(ltMinAge, ple$Ages[idpl])
            plev <- ple$ple[c(idpl[1] - 1, idpl)]
            plev <- plev / plev[1]
          } else {
            pleAges <- ple$Ages[idpl]
            plev <- ple$ple[idpl]
            plev <- plev / plev[1]
          }
        } else {
          pleAges <- ple$Ages
          plev <- ple$ple
        }
        lines(pleAges, plev, col= 'grey60', type = 's')
      }
      lines(xx + minAge, yy[1, ], lwd = lwdd, col = Palette[1], 
            lty = ltyy)
      if (icat == ncat) {
        legend('topright', c("Life table", "Estimated"), 
               col = c(1, Palette[1]), lwd = 2, bty = 'n', pch = c(19, NA))
      }
      
      # Mortality:
      xx <- x$x[cuts]
      yy <- x$mort[[icat]][, cuts]
      dxLt <- diff(lifeTab$Ages[1:2])
      # ======= BUG 2023-10-19 ========= #
      # ltMu <- -log(1 - lifeTab$qx)
      ltMu <- -log(1 - lifeTab$qx) / dxLt
      # ================================ #
      ltMu[which(is.infinite(ltMu))] <- NA
      xlimMort <- c(0, max(lifeTab$Ages))
      ylimMort <- c(0, max(c(ltMu), na.rm = TRUE))
      plot(xlimMort, ylimMort, col = NA, xlab = "Age", ylab = "Mortality")
      if (minAge > 0) lines(rep(minAge, 2), ylimMort, lty = 2, col = 'orange')
      if (!noCIs) {
        polygon(c(xx, rev(xx)) + minAge, c(yy[2, ], rev(yy[3, ])), 
                col = adjustcolor(Palette[1], alpha.f = 0.25),
                border = NA)
      }
      lines(xx + minAge, yy[1, ], col = Palette[1], lty = ltyy, lwd = lwdd)
      # ======= BUG 2023-10-19 ========= #
      lines(lifeTab$Ages[idAges] + dxLt / 2, ltMu[idAges], pch = 19, type = 'b',
            lty = 2)
      # lines(lifeTab$Ages[idAges], ltMu[idAges], pch = 19, type = 'b',
      #       lty = 2)
      # ================================ #
      
    }
  } else if (type == "fancy") {
    allThetaNames <- c("a0", "a1", "c", "b0", "b1", "b2")
    allThetaExpr  <- expression(italic(a[0]), italic(a[1]), italic(c), 
                                italic(b[0]), italic(b[1]), italic(b[2]))
    parNames <- substr(colnames(x$params), 1, 2)
    idTh <- which(is.element(substr(parNames, 1, 1), c("a", "b", "c")))
    thetaMat <- x$params[,idTh]
    model <- as.character(x$modelSpecs['model'])
    shape <- as.character(x$modelSpecs['shape'])
    
    if (model == "EX" | model == 1) {
      nTheta <- 1
      idAllTheta <- 4
    } else if (model %in% c("GO", "WE") | model %in% c(2, 3)) {
      nTheta <- 2
      idAllTheta <- c(4, 5)
    }  else {
      nTheta <- 3
      idAllTheta <- c(4, 5, 6)
    }
    if (shape == "Makeham" | shape == 2) {
      nTheta <- nTheta + 1
      idAllTheta <- c(3, idAllTheta)
    } else if(shape == "bathtub" | shape == 3) {
      nTheta <- nTheta + 3
      idAllTheta <- c(1:3, idAllTheta)
    }
    lows <- c(-Inf, 0, 0, -Inf, 0, 0)
    names(lows) <- allThetaNames
    
    # Mortality ylim:
    ylimm <- c(0, NA)
    for (icat in 1:lenCat) {
      ylimm[2] <- max(c(ylimm[2], x$mort[[icat]][, x$cuts[[icat]]]), 
                      na.rm = TRUE)
    }
    ylimm[2] <- ceiling(ylimm[2] * 10) / 10
    
    # Demo variable labels:
    demLab <- c(surv = expression(paste("Survival, ", italic(S(x)))),
                mort = expression(paste("Mortality, ", italic(mu(x)))))
    
    # Build plots:  
    layout(mat = matrix(data = c(rep(1:nTheta, each = 2), 
                                 rep(rep(c(nTheta + 1, nTheta + 2), 
                                         each = nTheta), 2)), 
                        nrow  = nTheta * 2, 
                        ncol  = 3), 
           widths  = rep(2, 3), 
           heights = rep(1, nTheta))
    par(mar=c(3, 3, 0.5, 0.5))
    for(i in 1:nTheta) {
      dez <- list()
      ylz <- rep(NA, lenCat)
      xlz <- matrix(0, lenCat, 2)
      for(j in 1:lenCat){
        idth <- paste(allThetaNames[idAllTheta[i]],
                      catNames[j], sep = ".")
        if (lenCat == 1) {
          idth <- allThetaNames[idAllTheta[i]]
        }
        dez[[catNames[j]]] <- density(thetaMat[, idth])
        ylz[j] <- max(dez[[catNames[j]]]$y)
        #xlz[j,] <- range(dez[[j]]$x)
        xlz[j, ] <- quantile(thetaMat[, idth], c(0.01, 0.99))
      }
      xr <- range(xlz)
      xl <- c(floor(max(c(lows[idAllTheta[i]], xr[1])) * 10) / 10, 
              ceiling(xr[2] * 10) / 10)
      xd <- ceiling(diff(xl) * 10) / 10
      plot(x = dez[[1]], xlab = "", ylab = "", xlim = xl, ylim = c(0, max(ylz)), 
           lwd = 3, axes = FALSE, main = "", col = NA)
      for(icat in 1:lenCat) {
        polygon(x = c(dez[[icat]]$x, dez[[icat]]$x[1]), 
                y = c(dez[[icat]]$y, dez[[icat]]$y[1]), 
                col = adjustcolor(Palette[icat], alpha.f = 0.25), 
                border = Palette[icat], lwd = 1.5)
      }
      axis(side = 1, at = seq(xl[1], xl[2], length = 5), 
           line = 0.5, labels = NA, tcl = 0.4)
      axis(side = 1, at = seq(xl[1], xl[2], length = 3), lwd = NA)
      mtext(text = allThetaExpr[idAllTheta[i]], side = 2, line = 0, 
            at = max(ylz) * 0.8, las = 2, cex = 1.25)
    }
    
    # Plot survival probability:
    par(mar = c(4, 7, 0.5, 0.5))
    xv <- lapply(1:lenCat, function(idcovs) x$x[x$cuts[[idcovs]]])
    mxv <- ceiling(max(unlist(xv)) / 5) * 5
    
    for (idem in c("mort", "surv")) {
      if (idem == "surv") {
        ylmx <- c(0, 1)
      } else {
        ylmx <- ylimm
      }
      plot(x = c(0, mxv), y = ylmx, col = NA, axes = FALSE, xlab = "", 
           ylab = "")
      for(icat in 1:lenCat) {
        xx <- x$x[x$cuts[[icat]]]
        yy <- x[[idem]][[icat]][, x$cuts[[icat]]]
        polygon(x = c(xx, rev(xx)), 
                y = c(yy[2, ], rev(yy[3, ])), 
                col = adjustcolor(Palette[icat], alpha.f = 0.25), 
                border = Palette[icat])
        lines(x = xx, y = yy[1, ], col = Palette[icat], lty = 3)
      }
      if (lenCat > 1 & idem == "surv") {
        legend(x = 'topright', legend = legNames, 
               pch = 22, pt.cex = 3, cex = 1.5, 
               pt.bg = adjustcolor(Palette, alpha.f = 0.25), col = Palette,
               bty = 'n')
      }
      
      axis(side = 2, at = seq(0, 1, 0.2), tcl = 0.4, las = 2, cex.axis = 1.2)
      mtext(text = demLab[idem], side = 2, line = 3.5, cex = 1.25)
      if (idem == "mort") {
        axis(side = 1, at = seq(0, mxv, 5), labels = NA, tcl = 0.4, line = 0.5)
      } else {
        axis(side = 1, at = seq(0, mxv, 5), tcl = 0.4, line = 0.5)
        mtext(text = expression(paste("Age ", italic(x), " (years)")), 
              side = 1, cex = 1.25, line = 3) 
      }
    }
  }
}

# Printing:
print.basta <- function(x, ...) {
  extraArgs <- list(...)
  if (length(extraArgs) > 0) {
    if (!is.element('digits', names(extraArgs))){
      digits <- 4
    } else {
      digits <- extraArgs$digits
    }
  } else {
    digits <- 4
  }
  if ("ModelSpecs" %in% names(x)) {
    x$modelSpecs <- x$ModelSpecs
  }
  cat("\nCall:\n")
  cat(paste("Model             \t\t: ", x$modelSpecs[1], "\n", sep = ""))
  cat(paste("Shape             \t\t: ", x$modelSpecs[2], "\n", sep = ""))
  cat(paste("Minimum age       \t\t: ", x$modelSpecs[3], "\n", sep = ""))
  cat(paste("Covars. structure \t\t: ", x$modelSpecs[4], "\n", sep = ""))
  cat(paste("Cat. covars.      \t\t: ", x$modelSpecs[5], "\n", sep = ""))
  cat(paste("Cont. covars.     \t\t: ", x$modelSpecs[6], "\n", 
            collapse = ""))
  
  cat("\nCoefficients:\n")
  print.default(x$coefficients, digits, ...)
  cat("\nConvergence:\n")
  cat(x$convmessage)
  if (is.na(x$DIC[1])){
    cat("\nDIC not calculated.")
  } else {
    cat(sprintf("\nDIC = %s\n", round(x$DIC["DIC"], 2)))
  }
}

# Summary:
summary.basta <- function(object, ...){
    extraArgs       <- list(...)
    if (length(extraArgs) > 0) {
      if (!is.element('digits', names(extraArgs))){
        digits <- 4
      } else {
        digits <- extraArgs$digits
      }
    } else {
      digits <- 4
    }
    if ("ModelSpecs" %in% names(object)) {
      object$modelSpecs <- object$ModelSpecs
    }
    if ("version" %in% names(object)) {
      cat(sprintf("\nOutput from BaSTA version %s\n", object$version))
    }
    cat("\nCall:\n")
    cat(paste("Model             \t\t: ", object$modelSpecs[1], "\n", sep = ""))
    cat(paste("Shape             \t\t: ", object$modelSpecs[2], "\n", sep = ""))
    cat(paste("Minimum age       \t\t: ", object$modelSpecs[3], "\n", sep = ""))
    cat(paste("Covars. structure \t\t: ", object$modelSpecs[4], "\n", sep = ""))
    cat(paste("Cat. covars.      \t\t: ", object$modelSpecs[5], "\n", sep = ""))
    cat(paste("Cont. covars.     \t\t: ", object$modelSpecs[6], "\n", 
              collapse = ""))
    
    cat("\nModel settings:\n")
    print(object$set)
    
    
    cat("\nMean Kullback-Leibler\ndiscrepancy calibration (KLDC):\n")
    if (object$K[1] != "Not calculated") {
      if ("qkl1" %in% names(object$K)) {
        meanKLcalib  <- t((object$K$qkl1 + object$K$qkl2) / 2)
      } else {
        meanKLcalib  <- (object$K$q12 + object$K$q21) / 2
      }
      print.default(meanKLcalib, digits = digits)
    } else {
      if (object$set['nsim'] == 1) {
        cat("KLDC was not calculated due to insufficient number\n",
            " of simulations to estimate convergence.\n")
      } else {
        cat("KLDC was not calculated due to lack of convergence,\n",
            "or because covariates were not included in the model.\n")
      }
    }
    
    
    cat("\nCoefficients:\n")
    print.default(object$coefficients, digits = digits)
    
    cat("\nConvergence:\n")
    if ("Convergence" %in% names(object)) {
      object$convergence <- object$Convergence
    }
    if (is.null(object$convergence)) {
      if (object$set['nsim'] == 1) {
        message("\nConvergence calculations require more than one run.",
            "\nTo estimate potential scale reduction run at least",
            " two simulations.\n")
      } else {
        cat("\nWarning: Convergence not reached for some parameters",
            " (i.e. 'PotScaleReduc' values larger than 1.1).",
            "\nThese estimates should not be used for inference.\n")
      }
    } else {
      if (all(object$convergence[, "Rhat"] < 1.1)) {
        cat("Appropriate convergence reached for all parameters.\n")
      } else {
        cat("\nWarning: Convergence not reached for some parameters",
            " (i.e. 'PotScaleReduc' values larger than 1.1).",
            "\nThese estimates should not be used for inference.\n")
      }
    } 
    cat("\nDIC:\n")
    if (!is.na(object$DIC[1])){
      cat(object$DIC["DIC"],"\n")
      if ("Convergence" %in% names(object)) {
        warning("Model fit in versions older than BaSTA 1.5 had a mistake in the",
                "\ncalculation of DIC values. In case you are interested in\n",
                "comparing the fit of different models, please try to run them\n",
                "with BaSTA versions 1.5 or above.",
                " We apologise for the inconvenience.", call. = FALSE)
      }
    } else {
      if (object$set['nsim'] == 1) {
        cat("DIC was not calculated due to insufficient number",
            "of simulations to estimate convergence.\n")
      } else {
        cat("DIC was not calculated due to lack of convergence.\n")
      }
    }
    ans <- c(list(coefficients = object$coef, DIC = object$DIC,
                  KullbackLeibler = object$KullbackLeibler, 
                  convergence = object$convergence,
                  modelSpecs = object$modelSpecs, settings = object$set))
    return(invisible(ans))
  }

# A.4) construct capture-recapture matrix:
# ---------------------------------------- #
CensusToCaptHist <- function(ID, d, dformat = "%Y", timeInt = "Y") {
  # Check data
  if(!inherits(ID, "character")) {
    ID <- as.character(ID)
  } 
  if (is.numeric(d)) {
    if (length(which(round(d) != d)) > 0) {
      stop("Please provide integer values or Date class ", 
           "values for argument 'd'.", call. = FALSE)
    } else {
      int <- d
    }
  } else if (is.character(d) | inherits(d, "Date")) {
    if (is.character(d)) {
      d <- as.Date(d, format = dformat)
      if (length(which(is.na(d)))) {
        stop("Wrong 'dformat' argument or wrong 'd' values.", 
             call. = FALSE)
      }
    }
    if (timeInt == "Y"){
      int <- as.numeric(format(d, format = "%Y"))
    } else if (timeInt == "M") {
      int <- as.numeric(format(d, format = "%m")) + 
        12 * (as.numeric(format(d, format = "%Y")) - 
                min(as.numeric(format(d, format = "%Y"))))
    } else if (timeInt == "D" | timeInt == "W") {
      jul <- julian(d, origin = min(d)) + 1
      if (timeInt == "W"){
        int <- ceiling(jul / 7)
      } else {
        int <- jul
      }
    }
  } else {
    stop("Wrong class for argument 'd'. Values\n", 
         "need to be of class 'integer', 'character' or 'Date'.", 
         call. = FALSE)
  }  
  
  # Construct capture-recapture matrix:
  dint <- min(int):max(int)
  ndint <- length(dint)
  uniID <- sort(unique(ID))
  n <- length(uniID)
  mat <- matrix(0, n, ndint, dimnames = list(NULL, dint))
  for (i in 1:ndint) {
    idt <- which(d == dint[i])
    idd <- which(uniID %in% ID[idt])
    mat[idd, i] <- 1
  }
  
  # Create data.frame with ID and Y:
  dmat <- data.frame(ID = uniID, mat)
  
  return(dmat)
}

# A.7) Multi-BaSTA functions:
# --------------------------- #
# Main multibasta:
multibasta <- function(object, dataType = "CMR", models, 
                       shapes, ...) {
  # Extract arguments:
  argList <- list(...)
  argNames <- names(argList)
  
  # Models available:
  mods <- rbind(data.frame(model = "EX", shape = "simple"),
                expand.grid(model = c("GO", "WE", "LO"), 
                            shape = c("simple", "Makeham", "bathtub")))
  
  # Abreviations for models and shapes:
  modsAbr <- rbind(data.frame(model = "Ex", shape = "Si"),
                   expand.grid(model = c("Go", "We", "Lo"), 
                               shape = c("Si", "Ma", "Bt")))
  
  # Alternative is models and shape are missing (i.e. run all):
  if (missing(models)) models <- c("EX", "GO", "WE", "LO")
  if (missing(shapes)) shapes <- c("simple", "Makeham", "bathtub")
  
  # Models request by user:
  idIncl <- which(mods$model %in% models & mods$shape %in% shapes)
  mods <- mods[idIncl, ]
  modsAbr <- modsAbr[idIncl, ]
  nmod <- nrow(mods)
  modNames <- paste(modsAbr$model, modsAbr$shape, sep = ".")

  # Define nsim:
  if (!"nsim" %in% argNames) {
    argList$nsim <- 2
  }
  
  # Logical whether DICs can be calculated:
  if (argList$nsim > 1) DIC <- TRUE else DIC <- FALSE
  
  # DIC matrix:
  if (DIC) {
    dics <- matrix(NA, nmod, 5, dimnames = 
                     list(modNames, c("D.ave", "D.mode", "pD", "k", "DIC")))
  }
  
  # list for models ran:
  runList <- list()
  
  # Run models:
  for (mod in 1:nmod) {
    updt <- sprintf("Run number %s, model: %s", mod, modNames[mod])
    cat(paste("\n", paste(rep("-", nchar(updt)), collapse = ""), sep = ""))
    cat(sprintf("\n%s\n", updt))
    cat(paste(paste(rep("-", nchar(updt)), collapse = ""), "\n", sep = ""))
    out <- basta(object, dataType = dataType, model = mods$model[mod], 
                 shape = mods$shape[mod], ...)
    runList[[modNames[mod]]] <- out
    if (DIC & !is.na(out$DIC[1])) {
      dics[mod, ] <- out$DIC
    }
  }
  
  # Fill up DIC matrix:
  if (DIC) {
    idNa <- which(is.na(dics[, "DIC"]))
    idnNa <- which(!is.na(dics[, "DIC"]))
    idModFit <- c(idnNa[sort.int(dics[idnNa, 5], index.return = TRUE)$ix], 
                  idNa)
    if (is.na(idModFit[1])) idModFit <- 0
    modFit <- data.frame(mods, dics)[idModFit, ]
    DICdiff <- modFit$DIC - modFit$DIC[1]
    Rank <- 1:nmod
    modFit <- data.frame(modFit, DICdiff, Rank)
  } else {
    modFit <- data.frame(DICs = "Not Calculated")
    idModFit <- NA
  }
  
  # Results list:
  results <- list(runs = runList[idModFit], DICs = modFit, models = mods)
  class(results) <- "multibasta"
  return(results)
}

# Printing:
print.multibasta <- function(x, ...) {
  argList <- list(...)
  if (is.null(argList$digits)) digits <- 3 else digits <- argList$digits
  cat("\nDICs:\n")
  print(x$DICs, digits = digits)
}

# Summary:
summary.multibasta <- function(object, ...) {
  argList <- list(...)
  if (is.null(argList$digits)) digits <- 3 else digits <- argList$digits
  if ("version" %in% names(object$runs[[1]])) {
    cat(sprintf("\nBaSTA version %s\n", object$runs[[1]]$version))
  }
  
  cat("\nCall:\n")
  cat(paste("Covars. structure \t\t: ", object$runs[[1]]$modelSpecs[3], "\n", 
            sep = ""))
  cat(paste("Minimum age       \t\t: ", object$runs[[1]]$modelSpecs[4], "\n", 
            sep = ""))
  cat(paste("Cat. covars.      \t\t: ", object$runs[[1]]$modelSpecs[5], "\n", 
            sep = ""))
  cat(paste("Cont. covars.     \t\t: ", object$runs[[1]]$modelSpecs[6], "\n", 
            collapse = ""))
  
  cat("\nModel settings:\n")
  print(object$runs[[1]]$set)
  
  cat("\nDICs:\n")
  print(object$DICs, digits = digits)
}

# Extract coefficients:
coef.multibasta <- function(object, showAll = FALSE, ...) {
  argList <- list(...)
  if (is.null(argList$digits)) digits <- 3 else digits <- argList$digits
  ans <- list()
  if (inherits(showAll,  "logical")) {
    if (showAll) nmod <- 1:length(object$runs) else nmod <- 1
  } else if (inherits(showAll, "numeric")) {
    if (length(showAll) == 1) {
      nmod <- 1:showAll
    } else {
      nmod <- showAll
    }
    nmod <- nmod[nmod <= length(object$runs)]
  } else {
    stop("\nArgument 'showAll' should be logical or an integer.\n", 
         call. = FALSE)
  }
  modNames <- names(object$runs)
  for (mod in nmod) {
    updt <- sprintf("DIC rank %s, model: %s", as.character(object$DICs$Rank)[mod], 
                    modNames[mod])
    cat(paste("\n", paste(rep("-", nchar(updt)), collapse = ""), sep = ""))
    cat(sprintf("\n%s\n", updt))
    cat(paste(paste(rep("-", nchar(updt)), collapse = ""), "\n", sep = ""))
    print.default(object$runs[[mod]]$coefficients, digits = digits)
    ans$coefficients[[modNames[mod]]] <- object$runs[[mod]]$coefficients
  }
  ans$DICs <- object$DICs
  return(invisible(ans))
}

# ============================================ #
# ==== B) FUNCTIONS FOR INTERNAL OBJECTS: ==== 
# ============================================ #
# Algorithm object function:
.CreateAlgObj <- function(dataType, model, shape, studyStart, studyEnd, 
                          minAge, covarsStruct, formulaMort, formulaRecap, 
                          recaptTrans, niter, burnin, thinning, updateJumps,
                          nsim, negSenescence) {
  return(list(dataType = dataType, model = model, shape = shape, 
              start = studyStart, end = studyEnd, minAge = minAge, 
              covStruc = covarsStruct, formulaMort = formulaMort, 
              formulaRecap = formulaRecap, recap = recaptTrans, 
              niter = niter, burnin = burnin, thinning = thinning, 
              updJump = updateJumps, nsim = nsim,
              negSenescence = negSenescence))
}

# Prepare data object:
.CreateDataObj <- function(object, algObj) {
  classDataObj <- c("bastacmr", "ageUpd")
  # Data Object for CMR data type:
  if (algObj$dataType == "CMR") {
    dataObj <- list()
    # Extract study year sequence and length:
    dataObj$study <- algObj$start:algObj$end
    dataObj$studyLen <- length(dataObj$study)
    
    # Number of observations:
    dataObj$n <- nrow(object)
    
    # Recapture matrix:
    Y <- as.matrix(object[, 1:dataObj$studyLen + 3])
    dataObj$Y <- Y
    dataObj$Y[Y > 1] <- 1
    colnames(dataObj$Y) <- dataObj$study
    
    # NOTE: addition of censTime for studies where individuals are censored
    #       before the end of the study (e.g. two studies of different 
    #       duration).
    # Find possible times of censoring before the study end:
    censTime <- rep(algObj$end, dataObj$n)
    for (tt in 1:dataObj$studyLen) {
      idCens <- which(Y[, tt] > 1)
      censTime[idCens] <- dataObj$study[tt]
    }
    dataObj$censTime <- censTime
    
    # Birth - death matrix:
    bd <- as.matrix(object[, 2:3])
    dataObj$bi <- bd[, 1]
    dataObj$di <- bd[, 2]
    bi0 <- which(dataObj$bi == 0 | is.na(dataObj$bi))
    if (length(bi0) > 0) {
      dataObj$idNoB <- bi0
      dataObj$updB <- TRUE
      dataObj$nUpdB <- length(bi0)
    } else {
      dataObj$updB <- FALSE
      dataObj$nUpdB <- 0
    }
    di0 <- which(dataObj$di == 0 | is.na(dataObj$di))
    if (length(di0) > 0) {
      dataObj$idNoD <- di0
      dataObj$updD <- TRUE
      dataObj$nUpdD <- length(di0)
    } else {
      dataObj$updD <- FALSE
      dataObj$nUpdD <- 0
    }
    
    if (!dataObj$updB & !dataObj$updD) {
      classDataObj[2] <- "noAgeUpd"
      dataObj$updA <- FALSE
    } else {
      dataObj$idNoA <- sort(unique(c(dataObj$idNoB, dataObj$idNoD)))
      dataObj$nUpdA <- length(dataObj$idNoA)
      dataObj$updA <- TRUE
      
      # NOTE: addition of min-max birth and death (2022-05-17):
      if ("Min.Birth" %in% colnames(object)) {
        dataObj$minBirth <- object[, "Min.Birth"]
        idMinB <- which(!is.na(dataObj$minBirth[dataObj$idNoB]))
        if (length(idMinB) > 0) {
          dataObj$idMinB <- dataObj$idNoB[idMinB]
          dataObj$updMinB <- TRUE
        } else {
          dataObj$idMinB <- NA
          dataObj$updMinB <- FALSE
        }
      } else {
        dataObj$minBirth <- rep(NA, dataObj$n)
        dataObj$idMinB <- NA
        dataObj$updMinB <- FALSE
      }
      if ("Max.Birth" %in% colnames(object)) {
        dataObj$maxBirth <- object[, "Max.Birth"]
        idMaxB <- which(!is.na(dataObj$maxBirth[dataObj$idNoB]))
        if (length(idMaxB) > 0) {
          dataObj$idMaxB <- dataObj$idNoB[idMaxB]
          dataObj$updMaxB <- TRUE
        } else {
          dataObj$idMaxB <- NA
          dataObj$updMaxB <- FALSE
        }
      } else {
        dataObj$maxBirth <- rep(NA, dataObj$n)
        dataObj$idMaxB <- NA
        dataObj$updMaxB <- FALSE
      }
      if ("Min.Death" %in% colnames(object)) {
        dataObj$minDeath <- object[, "Min.Death"]
        idMinD <- which(!is.na(dataObj$maxDeath[dataObj$idNoD]))
        if (length(idMinD) > 0) {
          dataObj$idMinD <- dataObj$idNoD[idMinD]
          dataObj$updMinD <- TRUE
        } else {
          dataObj$idMinD <- NA
          dataObj$updMinD <- FALSE
        }
      } else {
        dataObj$minDeath <- rep(NA, dataObj$n)
        dataObj$idMinD <- NA
        dataObj$updMinD <- FALSE
      }
      if ("Max.Death" %in% colnames(object)) {
        dataObj$maxDeath <- object[, "Max.Death"]
        idMaxD <- which(!is.na(dataObj$maxDeath[dataObj$idNoD]))
        if (length(idMaxD) > 0) {
          dataObj$idMaxD <- dataObj$idNoD[idMaxD]
          dataObj$updMaxD <- TRUE
        } else {
          dataObj$idMaxD <- NA
          dataObj$updMaxD <- FALSE
        }
        
      } else {
        dataObj$maxDeath <- rep(NA, dataObj$n)
        dataObj$idMaxD <- NA
        dataObj$updMaxD <- FALSE
      }
      
      # 4.1.2 Calculate first and last time observed 
      #       and total number of times observed:
      ytemp <- t(t(dataObj$Y) * dataObj$study)
      dataObj$lastObs <- c(apply(ytemp, 1, max))
      ytemp[ytemp == 0] <- 10000
      dataObj$firstObs <- c(apply(ytemp, 1, min))
      dataObj$firstObs[dataObj$firstObs == 10000] <- 0
      dataObj$oi <- dataObj$Y %*% rep(1, dataObj$studyLen)
      
      # 4.1.3 Define study duration:
      dataObj$Tm <- matrix(dataObj$study, dataObj$n, 
                           dataObj$studyLen, byrow = TRUE)
      fii <- dataObj$firstObs
      id1 <- which(dataObj$bi > 0 & dataObj$bi >= algObj$start)
      fii[id1] <- dataObj$bi[id1] + 1
      fii[dataObj$bi > 0 & dataObj$bi < algObj$start]  <- algObj$start
      lii <- dataObj$lastObs
      id2 <- which(dataObj$di > 0 & dataObj$di <= dataObj$censTime)
      lii[id2] <- dataObj$di[id2] - 1
      idCens <- which(dataObj$di > 0 & dataObj$di > dataObj$censTime)
      lii[idCens] <- dataObj$censTime[idCens]
      dataObj$obsMat <- .BuildAliveMatrix(fii, lii, dataObj)
      dataObj$obsMat[lii == 0 | fii == 0, ] <- 0
      classDataObj[2] <- "ageUpd"
    }
    dataObj$Dx <- 1 
    classDataObj[1] <- "bastacmr"
  }
  
  # Data Object for Census data type:
  else {
    n <- nrow(object)
    # Calculate Julian times:
    # bi <- round(as.numeric(as.Date(object$Birth.Date, 
    #                                format = "%Y-%m-%d")) / 
    #               365.25, 2) + 1970
    # bil <- round(as.numeric(as.Date(object$Min.Birth.Date, 
    #                                 format = "%Y-%m-%d")) /
    #                365.25, 2) + 1970
    # biu <- round(as.numeric(as.Date(object$Max.Birth.Date, 
    #                                 format = "%Y-%m-%d")) /
    #                365.25, 2) + 1970
    # firstObs <- round(as.numeric(as.Date(object$Entry.Date, 
    #                                   format = "%Y-%m-%d")) /
    #                  365.25, 2) + 1970
    # lastObs <- round(as.numeric(as.Date(object$Depart.Date, 
    #                                    format = "%Y-%m-%d")) /
    #                   365.25, 2) + 1970

    bi <- as.numeric(as.Date(object$Birth.Date, 
                             format = "%Y-%m-%d")) / 365.25 + 1970
    bil <- as.numeric(as.Date(object$Min.Birth.Date, 
                              format = "%Y-%m-%d")) / 365.25 + 1970
    biu <- as.numeric(as.Date(object$Max.Birth.Date, 
                              format = "%Y-%m-%d")) / 365.25 + 1970
    firstObs <- as.numeric(as.Date(object$Entry.Date, 
                                   format = "%Y-%m-%d")) / 365.25 + 1970
    lastObs <- as.numeric(as.Date(object$Depart.Date, 
                                  format = "%Y-%m-%d")) / 365.25 + 1970
    
    # Entry and departure types:
    entryType <- as.character(object$Entry.Type)
    departType <- as.character(object$Depart.Type)
    
    # Censored individuals:
    idCens <- which(departType %in% c("O", "C"))
    nCens <- length(idCens)
    
    # Birth times to be estimated:
    idNoBirth <- which(bi != bil | bi != biu)
    nNoBirth <- length(idNoBirth)
    if (nNoBirth == 0) {
      updB <- FALSE
      classDataObj[2] <- "noAgeUpd"
    } else {
      updB <- TRUE
      classDataObj[2] <- "ageUpd"
    }

    # Create data object:
    dataObj <- list(bi = bi, bil = bil, biu = biu, firstObs = firstObs, 
                    lastObs = lastObs, idCens = idCens,
                    idNoB = idNoBirth, nUpdB = nNoBirth, updB = updB, 
                    updD = FALSE, idNoD = NA, nUpdD = 0,
                    idNoA = idNoBirth, nUpdA = nNoBirth, updA = updB, n = n, 
                    nCens = nCens, studyLen = NA)
    classDataObj[1] <- "bastacensus"
  }
  class(dataObj) <- classDataObj

  # Output:
  return(dataObj)
}

# Create initial age object:
.CreateAgeObj <- function(dataObj, algObj) {
  ageObj <- list()
  bi <- dataObj$bi
  if (dataObj$updB & inherits(dataObj, "bastacmr")) {
    bi[dataObj$idNoB] <- dataObj$firstObs[dataObj$idNoB] - 
      sample(6:1, size = dataObj$nUpdB, replace = TRUE)
  }
  if (inherits(dataObj, "bastacmr")) {
    di <- dataObj$di
    if (dataObj$updD) {
      di[dataObj$idNoD] <- apply(cbind(bi[dataObj$idNoD], 
                                       dataObj$lastObs[dataObj$idNoD]),
                                 1, max) + 
        sample(6:1, size = dataObj$nUpdD, replace = TRUE)
    }
    age <- di - bi
    ageTr <- algObj$start - bi
    ageTr[ageTr < 0] <- 0
  } else {
    di <- dataObj$lastObs
    age <- dataObj$lastObs - bi
    ageTr <- dataObj$firstObs - bi
    ageTr[ageTr < 0] <- 0
  }
  
  # indicator for uncensored:
  indUncens <- rep(1, dataObj$n)
  if (inherits(dataObj, "bastacensus")) {
    indUncens <- rep(1, dataObj$n)
    indUncens[dataObj$idCens] <- 0
  }
  
  # Recalculate ages based on minAge:
  # ---------------------------------
  if (algObj$minAge > 0) {
    # Ages and indicator after min age:
    ageAftMa <- age - algObj$minAge
    ageAftMa[ageAftMa < 0] <- 0
    
    # Ages and indicator before min age:
    ageBefMa <- age
    indBefMa <- rep(0, dataObj$n)
    indBefMa[ageBefMa < algObj$minAge] <- 1
    ageBefMa[age >= algObj$minAge] <- algObj$minAge
    
    # Ages at truncation and indicator after min age:
    ageTrAftMa <- ageTr - algObj$minAge
    ageTrAftMa[ageTrAftMa < 0] <- 0
    
    # Ages at truncation and indicator before min age:
    ageTrBefMa <- ageTr
    ageTrBefMa[ageTr >= algObj$minAge] <- algObj$minAge
  } else {
    ageBefMa <- ageTrBefMa <- indBefMa <- rep(0, dataObj$n)
    ageAftMa <- age
    ageTrAftMa <- ageTr
  }
  
  
  # Create alive matrix for bastacmr:
  # ---------------------------------
  if (inherits(dataObj, "bastacmr")) {
    firstObs <- c(apply(cbind(algObj$start, bi + 1), 1, max))
    lastObs <- c(apply(cbind(algObj$end, dataObj$censTime, di), 1, min))
    alive <- .BuildAliveMatrix(firstObs, lastObs, dataObj)
  } else {
    alive <- NA
  }
  
  # Fill-in matrices for age object:
  # --------------------------------
  ageObj$ages <- data.frame(birth = bi, death = di, age = age, ageTr = ageTr,
                            ageAft = ageAftMa, ageBef = ageBefMa, 
                            truAft = ageTrAftMa, truBef = ageTrBefMa)
  ageObj$inds <- data.frame(ageBef = indBefMa, uncens = indUncens)
  ageObj$alive <- alive
  
  # Assign class to age object:
  # ---------------------------
  minAgeClass <- ifelse(algObj$minAge > 0, "minAge", "noMinAge")
  if (inherits(dataObj, "bastacmr")) {
    class(ageObj) <- c("agecmr", minAgeClass)
  } else {
    class(ageObj) <- c("agecensus", minAgeClass)
  }
  return(ageObj)
}

# Build a matrix with 1 when alive, 0 otherwise:
.BuildAliveMatrix <- function(f, l, dataObj) {
  Fm <- dataObj$Tm - f
  Fm[Fm >= 0] <- 1
  Fm[Fm < 0] <- 0
  Lm <- dataObj$Tm - l
  Lm[Lm <= 0] <- -1
  Lm[Lm > 0] <- 0
  return(Fm * (-Lm))
}

# Create covariates object:
.CreateCovObj <- function(object, dataObj, algObj) {
  covObj <- list()
  # ------------------------- #
  # Covariates for mortality:
  # ------------------------- #
  covClass <- c("noCov", "noCovType", "noCovRecap")
  # Find if mortality covariates are not needed:
  # -------------------------------------------- #
  if (is.null(algObj$formulaMort)) {
    covObj$covs <- NULL
  } else {
    # When there should be covariates for mortality:
    # ---------------------------------------------- #
    # Construct covariate matrix:
    covMat <- .MakeCovMat(algObj$formulaMort, data = object)
    
    # Find covariate types (categorical vs continuous):
    covType <- .FindCovType(covMat)
    
    # Separate covariates by their type and by covarsStruct:
    if (algObj$covStruc == "fused") {
      covClass[1] <- "fused"
      if (!is.null(covType$cat)) {
        covObj$inMort <- covMat[, covType$cat]
        covObj$imLen <- ncol(covObj$inMort)
      } else {
        covClass[1] <- "propHaz"
      }
      if (!is.null(covType$cont)) {
        covObj$propHaz <- matrix(covMat[, c(covType$int, covType$cont)], 
                                 ncol = length(c(covType$int, covType$cont)),
                                 dimnames = list(NULL, c(names(covType$int), 
                                                         names(covType$cont))))
        covObj$phLen <- ncol(covObj$propHaz)
      } else {
        covClass[1] <- "inMort"
      }
    } else if (algObj$covStruc == "all.in.mort") {
      if (is.null(covType$int) & is.null(covType$cat)) {
        covObj$inMort <- cbind(1, covMat)
        colnames(covObj$inMort) <- c("Intercept", colnames(covMat))
      } else {
        covObj$inMort <- covMat
      }
      covObj$imLen <- ncol(covObj$inMort)
      covClass[1] <- "inMort"
    } else {
      if (!is.null(covType$int)) {
        covObj$propHaz <- 
          matrix(covMat[, -covType$int], dataObj$n, ncol(covMat) -1, 
                 dimnames = list(NULL, colnames(covMat)[-covType$int]))
      } else if (!is.null(covType$cat)) {
        covObj$propHaz <- 
          matrix(covMat[, -covType$cat[1]], dataObj$n, ncol(covMat) -1, 
                 dimnames = list(NULL, colnames(covMat)[-covType$cat[1]]))
      } else {
        covObj$propHaz <- covMat
      }
      covObj$phLen <- ncol(covObj$propHaz)
      covClass[1] <- "propHaz"
    }
    if (!is.null(covType$cat) & !is.null(covType$cont)) {
      covClass[2] <- "bothCov"
      covObj$cat <- covType$cat
      covObj$cont <- covType$cont
    } else if (!is.null(covType$cat)) {
      covClass[2] <- "cateCov"
      covObj$cat <- covType$cat
    } else if (!is.null(covType$cont)) {
      covClass[2] <- "contCov"
      covObj$cont <- covType$cont
    }
    
    # General covariate names:
    genCovs <- strsplit(gsub(" - 1", "", as.character(algObj$formulaMort)[2]),
                        "[[:blank:]]{1}[[:punct:]]{1}[[:blank:]]{1}")[[1]]
    covObj$genCovs <- genCovs
  }
  
  # ------------------------------------- #
  # Covariates for recapture probability:
  # ------------------------------------- #
  
  # ------------------------ #
  # Return covariate object:
  # ------------------------ #
  class(covObj) <- covClass
  return(covObj)
}

# Find covariate type (i.e. continuous vs categorical):
.FindCovType <- function(covMat) {
  # This functions finds and returns if an intercecpt was included 
  # and which covariates are categorical or continuous.
  if (!is.null(covMat)) {
    lu <- apply(covMat, 2, function(x) length(unique(x)))
    ru <- apply(covMat, 2, range)
    idcat <- which(lu == 2 & apply(ru, 2, sum) == 1)
    if (length(idcat) == 0) {
      idcat <- NULL
    }
    idint <- which(lu == 1)
    if (length(idint) == 0) {
      idint <- NULL
    }
    idcon <- which(lu > 2)
    if (length(idcon) == 0) {
      idcon <- NULL
    }
  }
  else {
    idcat <- NULL
    idint <- NULL
    idcon <- NULL
  }
  return(list(int = idint, cat = idcat, cont = idcon))
}

# Construct a design matrix (i.e. covariate matrix) following a formula:
.MakeCovMat <- function(covform, data) {
  covs <- model.matrix(covform, data = data)
  return(covs)
}

# ------------------------------- #
# Functions to manage parameters:
# ------------------------------- #
# Extract parameters specified by the user:
.CreateUserPar <- function(covObj, argList) {
  userPars <- list()
  genParName <- c("theta", "gamma")
  parTypes <- c("Start", "Jumps", "PriorMean", "PriorSd")
  parTypesList <- c("start", "jump", "priorMean", "priorSd")
  for (genPp in 1:2) {
    if (all(genPp == 2 & inherits(covObj, c("fused", "propHaz"))) |
        genPp == 1) {
      userPars[[genParName[genPp]]] <- list()
      for (pp in 1:4) {
        usrPar <- sprintf("%s%s", genParName[genPp], parTypes[pp])
        if (usrPar %in% names(argList)) {
          userPars[[genParName[genPp]]][[parTypesList[pp]]] <- 
            argList[[usrPar]]
        } else {
          userPars[[genParName[genPp]]][[parTypesList[pp]]] <- NULL 
        }
      }
    }    
  }
  return(userPars)
}

# Set the default mortality parameters:
.SetDefaultTheta  <- function(algObj) {
  if (algObj$model == "EX") {
    nTh <- 1
    startTh <- 0.2 
    jumpTh <- 0.1
    priorMean <- 0.06
    priorSd <- 1
    nameTh <- "b0"
    lowTh <- 0
    jitter <- 0.5
  } else if (algObj$model == "GO") {
    nTh <- 2 
    startTh <- c(-2, 0.01) 
    jumpTh <- c(0.1, 0.1)
    priorMean <- c(-3, 0.01)
    priorSd <- c(5, 1) # 2023-02-18
    nameTh <- c("b0", "b1")
    lowTh <- c(-Inf, 0)
    if (algObj$negSenescence) lowTh[2] <- -Inf
    jitter <- c(0.5, 0.2) 
    if (algObj$shape == "bathtub") {
      lowTh <- c(-Inf, 0)
    }
  } else if (algObj$model == "WE") {
    nTh <- 2
    startTh <- c(1.5, 0.2) 
    jumpTh <- c(.01, 0.1)
    priorMean <- c(1.5, .05)
    priorSd <- c(1, 1)
    nameTh <- c("b0", "b1")
    lowTh <- c(0, 0)
    jitter <- c(0.5, 0.2) 
  } else if (algObj$model == "LO") {
    nTh <- 3 
    startTh <- c(-2, 0.01, 1e-04) 
    jumpTh <- c(0.1, 0.1, 0.1) 
    priorMean <- c(-3, 0.01, 1e-10)
    priorSd <- c(1, 1, 1)
    nameTh <- c("b0", "b1", "b2")
    lowTh <- c(-Inf, 0, 0)
    jitter <- c(0.5, 0.2, 0.5) 
  }
  if (algObj$shape == "Makeham") {
    nTh <- nTh + 1 
    startTh <- c(0, startTh) 
    jumpTh <- c(0.1, jumpTh) 
    priorMean <- c(0, priorMean)
    priorSd <- c(1, priorSd)
    nameTh <- c("c", nameTh)
    lowTh <- c(0, lowTh)
    jitter <- c(0.25, jitter) 
  } else if (algObj$shape == "bathtub") {
    nTh <- nTh + 3 
    startTh <- c(-0.1, 0.6, 0, startTh)
    jumpTh <- c(0.1, 0.1, 0.1, jumpTh) 
    priorMean <- c(-2, 0.01, 0, priorMean)
    priorSd <- c(1, 5, 1, priorSd) # 2023-02-18
    nameTh <- c("a0", "a1", "c", nameTh)
    lowTh <- c(-Inf, 0, 0, lowTh)
    jitter <- c(0.5, 0.2, 0.2, jitter) 
  }
  defaultTheta  <- list(length = nTh, start = startTh, jump = jumpTh, 
                        priorMean = priorMean, priorSd = priorSd, name = nameTh, 
                        lower = lowTh, jitter = jitter)
  attr(defaultTheta, "model") = algObj$model
  attr(defaultTheta, "shape") = algObj$shape
  return(defaultTheta)
}

# Build the parameters object:
.CreateFullParObj <- function(covObj, defTheta, algObj, 
                             userPars, dataObj) {
  fullParObj <- list()
  fullParObj$theta <- list()
  statName <- c("start", "priorMean", "priorSd", "jump", "lower", "jitter")
  nstat <- length(statName)
  for (st in 1:nstat) {
    if (inherits(covObj, c("inMort", "fused"))) {
      if (is.null(userPars$theta[[statName[st]]])) {
        thetaMat <- matrix(defTheta[[statName[st]]], covObj$imLen, 
                           defTheta$length, byrow = TRUE, 
                           dimnames = list(colnames(covObj$inMort), 
                                           defTheta$name))
        if (st %in% c(1, 2)) {
          if (inherits(covObj, "inMort") & 
              inherits(covObj, c("contCov", "bothCov"))) {
            thetaMat[names(covObj$cont), ] <- 0
          }
        }
        fullParObj$theta[[statName[[st]]]] <- thetaMat
      } else {
        if (is.element(length(userPars$theta[[statName[st]]]), 
                       c(defTheta$length, defTheta$length * covObj$imLen))) {
          if (length(userPars$theta[[statName[st]]]) == defTheta$length) {
            fullParObj$theta[[statName[[st]]]] <- 
              matrix(userPars$theta[[statName[st]]], covObj$imLen, 
                     defTheta$length, byrow = TRUE, 
                     dimnames = list(colnames(covObj$inMort), defTheta$name))
          } else {
            fullParObj$theta[[statName[[st]]]] <- userPars$theta[[statName[st]]]
            dimnames(fullParObj$theta[[statName[[st]]]]) <- 
              list(colnames(covObj$inMort), defTheta$name)
          }
        } else {
          stop(paste("\nDimensions of theta ", statName[st], 
                     " matrix are incorrect.\n",
                     "Provide a single vector of length ", defTheta$length,
                     "\nor a matrix of dimensions ", covObj$imLen ," times ", 
                     defTheta$length, 
                     ".\n(i.e. number of covariates times number", 
                     " of\n parameters for model ", 
                     algObj$model," with ", algObj$shape, " shape).", 
                     sep = ""), call. = FALSE)
        }
      }
      allstatName <- paste(rep(defTheta$name, 
                               each = ncol(covObj$inMort)), 
                           rep(colnames(covObj$inMort), defTheta$len), 
                           sep = ".")
    } else {
      if (is.null(userPars$theta[[statName[st]]])) {
        fullParObj$theta[[statName[[st]]]] <- 
          matrix(defTheta[[statName[st]]], 1, defTheta$length, 
                 dimnames = list(NULL, defTheta$name))
      } else {
        if (length(userPars$theta[[statName[st]]]) == defTheta$length) {
          fullParObj$theta[[statName[[st]]]] <- 
            matrix(userPars$theta[[statName[st]]], 1, defTheta$length, 
                   dimnames = list(NULL, defTheta$name))
        } else {
          stop(paste("\nLength of theta ", statName[st], " is incorrect.\n",
                     "Provide a single vector of length ", defTheta$length,
                     ".\n(i.e. number of parameters for model ", 
                     algObj$model," with ", algObj$shape, " shape).", 
                     sep = ""), call. = FALSE)
        }
      }
      allstatName <- defTheta$name
    }
  }
  fullParObj$theta$lower <- t(t(fullParObj$theta$start) * 0 + defTheta$lower)
  fullParObj$theta$len <- length(fullParObj$theta$start)
  fullParObj$theta$names <- allstatName
  if (inherits(covObj, c("propHaz", "fused"))) {
    fullParObj$gamma <- list()
    for (st in 1:nstat) {
      if (is.null(userPars$gamma[[statName[st]]])) {
        fullParObj$gamma[[statName[st]]] <- 
          rep(c(0.01, 0, 1, 0.1, -Inf, 0.2)[st], covObj$phLen)
        names(fullParObj$gamma[[statName[st]]]) <- colnames(covObj$propHaz)
      } else {
        if (length(userPars$gamma[[statName[st]]]) == covObj$phLen) {
          fullParObj$gamma[[statName[st]]] <- userPars$gamma[[statName[st]]]
          names(fullParObj$gamma[[statName[st]]]) <- colnames(covObj$propHaz)
        } else {
          stop(paste("\nLength of gamma parameters is incorrect.\n",
                     "Provide a single vector of length ", covObj$phLen,
                     ".\n(i.e. number of proportional hazards covariates).", 
                     sep = ""), call. = FALSE)
        }
      }
    }
    fullParObj$gamma$len <- length(fullParObj$gamma$start)
    allstatName <- c(allstatName, paste("gamma", colnames(covObj$propHaz), 
                                        sep = "."))
    fullParObj$gamma$names <- paste("gamma", colnames(covObj$propHaz), 
                                    sep = ".")
  }
  if (inherits(covObj, c("inMort", "noCov"))) {
    Classes <- "theta"
  } else {
    Classes <- "theGam"
  }
  # Minimum age's lambda:
  if (algObj$minAge > 0) {
    fullParObj$lambda <- list(start = 0.01, priorMean = 0.01, priorSd = 1, 
                              jump = 0.01, lower = 0, jitter = 0.1,  len = 1)
    Classes <- c(Classes, "lambda")
  } else {
    Classes <- c(Classes, "noLambda")
  }
  
  # Detection probability:
  if (inherits(dataObj,  "bastacmr")) {
    if (inherits(dataObj, "ageUpd")) {
      fullParObj$pi <- list()
      study <- algObj$start:algObj$end
      if (length(algObj$recap) == length(study)) {
        if (all(algObj$recap %in% study)) {
          idpi <- rep(1, length(study))
          for (ipi in 2:length(idpi)) {
            if (algObj$recap[ipi] == algObj$recap[ipi - 1]) {
              idpi[ipi] <- idpi[ipi - 1]
            } else if (algObj$recap[ipi] %in% algObj$recap[1:(ipi - 1)]) {
              idpi[ipi] <- idpi[which(algObj$recap[1:(ipi - 1)] ==
                                        algObj$recap[ipi])[1]]
            } else {
              idpi[ipi] <- max(idpi[1:(ipi - 1)] + 1)
            }
          }
        } else if (all(algObj$recap %in% 1:length(study))) {
          idpi <- algObj$recap
        }
        namespi <- unique(algObj$recap)
      } else {
        idpi <- findInterval(study, algObj$recap)
        namespi <- algObj$recap
      }
      names(idpi) <- study
      npi <- length(unique(idpi))
      fullParObj$pi$start <- rep(0.5, npi)
      names(fullParObj$pi$start) <- namespi
      fullParObj$pi$idpi <- idpi
      fullParObj$pi$n <- npi
      fullParObj$pi$prior2 <- 1
      fullParObj$pi$Prior1 <- tapply(1 + t(t(dataObj$Y) %*% rep(1, dataObj$n)),
                                     idpi, sum)
      fullParObj$pi$len <- length(fullParObj$pi$start)
      if (length(namespi) > 1) {
        piNames <- paste("pi", namespi, sep = ".")
      } else {
        piNames <- "pi"
      }
      fullParObj$pi$names <- piNames
      fullParObj$allNames <- c(allstatName, piNames)
      Classes <- c(Classes, "pi", "noEta")
    } else {
      fullParObj$allNames <- allstatName
      Classes <- c(Classes, "noPi", "noEta")
    }
  } else {
    Classes <- c(Classes, "noPi", "noEta")
  }
  fullParObj$class <- Classes
  class(fullParObj) <- Classes
  return(fullParObj)
}

# Define initial set of parameters:
.CreateParObj <- function(fullParObj) {
  parObj <- list()
  parObj$theta <- fullParObj$theta$start
  if (inherits(fullParObj,  "theGam")) {
    parObj$gamma <- fullParObj$gamma$start
  } else {
    parObj$gamma <- 0
  }
  if (inherits(fullParObj,  "lambda")) {
    parObj$lambda <- fullParObj$lambda$start
  } else {
    parObj$lambda <- 0
  }
  if (inherits(fullParObj,  c("pi", "piEta"))) {
    parObj$pi <- fullParObj$pi$start
  }
  class(parObj) <- class(fullParObj)
  return(parObj)
}

# Create and manage design matrices:
.CalcCovPars <- function(parObj, parCov, covObj, dataObj, type = "theta") {
  if (is.null(parCov)) {
    parCovObjn <- list()
  } else {
    parCovObjn <- parCov
  }
  if (type %in% c("theta", "both")) {
    if (inherits(covObj, c("fused", "inMort"))) {
      parCovObjn$theta <- covObj$inMort %*% parObj$theta
    } else {
      parCovObjn$theta <- matrix(1, nrow = dataObj$n) %*% parObj$theta
    }
  }
  if (type %in% c("gamma", "both")) {
    if (inherits(parObj, "theGam")) {
      parCovObjn$gamma <- covObj$propHaz %*% parObj$gamma
    } else {
      parCovObjn$gamma <- rep(0, dataObj$n)
    }
  }
  return(parCovObjn)
}

# =================================== #
# ==== C) STATISTICAL FUNCTIONS: ====
# =================================== #
# --------------------------------- #
# Truncated distribution functions:
# --------------------------------- #
# Truncated normal:
.rtnorm <- function(n, mean, sd, lower = -Inf, upper = Inf) {
  Flow <- pnorm(lower, mean, sd)
  Fup <- pnorm(upper, mean, sd)
  ru <- runif(n, Flow, Fup)
  rx <- qnorm(ru, mean, sd)
  return(rx)
}

.dtnorm <- function(x, mean, sd, lower = -Inf, upper = Inf, log = FALSE) {
  Flow <- pnorm(lower, mean, sd)
  Fup <- pnorm(upper, mean, sd)
  densx <- dnorm(x, mean, sd) / (Fup - Flow)
  if (log) densx <- log(densx)
  return(densx)
}

.ptnorm <- function(q, mean, sd, lower = -Inf, upper = Inf, log = FALSE) {
  p <- (pnorm(q, mean, sd) - pnorm(lower, mean, sd)) / 
    (pnorm(upper, mean, sd) - pnorm(lower, mean, sd))
  if (log) {
    p <- log(p)
  }
  return(p)
}

.qtnorm <- function (p, mean = 0, sd = 1, lower = -Inf, upper = Inf) {
  p2 <- (p) * (pnorm(upper, mean, sd) - pnorm(lower, mean, sd)) + 
    pnorm(lower, mean, sd)
  q <- qnorm(p2, mean, sd)
  return(q)
}

# ---------------------------------- #
# Basic survival analysis functions:
# ---------------------------------- #
# a) General mortality functions:
.DefineMortMatrix <- function(algObj) {
  if (algObj$model == "EX") {
    mortfun <- function(theta, x) {
      if (length(theta) == length(x)) {
        mux <- theta
      } else {
        mux <- rep(theta, length(x))
      }
      return(mux)
    }
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        exp(theta[ ,"b0"] + theta[, "b1"] * x)
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        theta[, "c"] + exp(theta[, "b0"] + theta[, "b1"] * x)
      }
    } else {
      mortfun <- function(theta, x) {
        exp(theta[, "a0"] - theta[, "a1"] * x) + theta[, "c"] + 
          exp(theta[, "b0"] + theta[, "b1"] * x)
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        mu <- theta[, "b0"] * theta[, "b1"]^theta[, "b0"] * 
          x^(theta[, "b0"] - 1)
        id0 <- which(x == 0)
        if (length(id0) > 0) mu[id0] <- min(mu[-id0])
        return(mu)
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        mu <- theta[, "c"] + theta[, "b0"] * theta[, "b1"]^theta[, "b0"] * 
          x^(theta[, "b0"] - 1)
        id0 <- which(x == 0)
        if (length(id0) > 0) mu[id0] <- theta[id0, "c"]
        return(mu)      
      }
    } else {
      mortfun <- function(theta, x) {
        mu <- exp(theta[, "a0"] - theta[, "a1"] * x) + theta[, "c"] + 
          theta[, "b0"] * theta[, "b1"]^theta[, "b0"] * 
          x^(theta[, "b0"] - 1)
        id0 <- which(x == 0)
        if (length(id0) > 0) {
          mu[id0] <- exp(theta[id0, "a0"]) + theta[id0, "c"]
        }
        return(mu)
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        exp(theta[, "b0"] + theta[, "b1"] * x) / 
          (1 + theta[, "b2"] * exp(theta[, "b0"]) / 
             theta[, "b1"] * (exp(theta[, "b1"] * x) - 1))
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        theta[, "c"] + exp(theta[, "b0"] + theta[, "b1"] * x) / 
          (1 + theta[, "b2"] * exp(theta[, "b0"]) / 
             theta[, "b1"] * (exp(theta[, "b1"] * x) - 1))
      }
    } else {
      mortfun <- function(theta, x) {
        exp(theta[, "a0"] - theta[, "a1"] * x) + theta[, "c"] + 
          exp(theta[, "b0"] + theta[, "b1"] * x) / 
          (1 + theta[, "b2"] * exp(theta[, "b0"]) / 
             theta[, "b1"] * (exp(theta[, "b1"] * x) - 1))
      }
    }
  }
  return(mortfun)
}

.DefineMortNumeric <- function(algObj) {
  if (algObj$model == "EX") {
    mortfun <- function(theta, x) {
      if (length(theta) == length(x)) {
        mux <- theta
      } else {
        mux <- rep(theta, length(x))
      }
      return(mux)
    }
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        exp(theta["b0"] + theta["b1"] * x)
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        theta["c"] + exp(theta["b0"] + theta["b1"] * x)
      }
    } else {
      mortfun <- function(theta, x) {
        exp(theta["a0"] - theta["a1"] * x) + theta["c"] + 
          exp(theta["b0"] + theta["b1"] * x)
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        mu <- theta["b0"] * theta["b1"]^theta["b0"] * 
          (x + 0.003)^(theta["b0"] - 1)
        return(mu)
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        mu <- theta["c"] + theta["b0"] * theta["b1"]^theta["b0"] * 
          (x + 0.003)^(theta["b0"] - 1)
        return(mu)
      }
    } else {
      mortfun <- function(theta, x) {
        mu <- exp(theta["a0"] - theta["a1"] * x) + theta["c"] + 
          theta["b0"] * theta["b1"]^theta["b0"] * 
          (x + 0.003)^(theta["b0"] - 1)
        id0 <- which(x == 0)
        if (length(id0) > 0) {
          mu[id0] <- exp(theta["a0"]) + theta["c"]
        }
        return(mu)
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      mortfun <- function(theta, x) {
        exp(theta["b0"] + theta["b1"] * x) / 
          (1 + theta["b2"] * exp(theta["b0"]) / 
             theta["b1"] * (exp(theta["b1"] * x) - 1))
      }
    } else if (algObj$shape == "Makeham") {
      mortfun <- function(theta, x) {
        theta["c"] + exp(theta["b0"] + theta["b1"] * x) / 
          (1 + theta["b2"] * exp(theta["b0"]) / 
             theta["b1"] * (exp(theta["b1"] * x) - 1))
      }
    } else {
      mortfun <- function(theta, x) {
        exp(theta["a0"] - theta["a1"] * x) + theta["c"] + 
          exp(theta["b0"] + theta["b1"] * x) / 
          (1 + theta["b2"] * exp(theta["b0"]) / 
             theta["b1"] * (exp(theta["b1"] * x) - 1))
      }
    }
  }
  return(mortfun)
}

# b) General cumulative hazards functions:
.DefineCumHazMatrix <- function(algObj) {
  if (algObj$model == "EX") {
    cumhazfun <- function(theta, x) c(theta) * x
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        exp(theta[, "b0"]) / theta[, "b1"] * 
          (exp(theta[, "b1"] * x) - 1)
      }
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta[, "c"] * x + exp(theta[, "b0"]) / theta[, "b1"] * 
          (exp(theta[, "b1"] * x) - 1)
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta[, "a0"]) / theta[, "a1"] * (1 - exp(-theta[, "a1"] * x)) + 
          theta[, "c"] * x + exp(theta[, "b0"]) / theta[, "b1"] * 
          (exp(theta[, "b1"] * x) - 1)
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        (theta[, "b1"] * x)^theta[, "b0"]
      }      
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta[, "c"] * x + (theta[, "b1"] * x)^theta[, "b0"]
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta[, "a0"]) / theta[, "a1"] * (1 - exp(-theta[, "a1"] * x)) +
          theta[, "c"] * x + (theta[, "b1"] * x)^theta[, "b0"]
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        log(1 + theta[, "b2"] * exp(theta[, "b0"]) / theta[, "b1"] * 
              (exp(theta[, "b1"] * x) - 1)) * (1 / theta[, "b2"])
      }
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta[, "c"] * x + log(1 + theta[, "b2"] * exp(theta[, "b0"]) / 
                                 theta[, "b1"] * 
                                 (exp(theta[, "b1"] * x) - 1)) * 
          (1 / theta[, "b2"])
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta[, "a0"]) / theta[, "a1"] * (1 - exp(-theta[, "a1"] * x)) +
          theta[, "c"] * x + log(1 + theta[, "b2"] * 
                                   exp(theta[, "b0"]) / theta[, "b1"] * 
                                   (exp(theta[, "b1"] * x) - 1)) *
          (1 / theta[, "b2"])
      }
    }
  }
  return(cumhazfun)
}

.DefineCumHazNumeric <- function(algObj) {
  if (algObj$model == "EX") {
    cumhazfun <- function(theta, x) c(theta) * x
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        exp(theta["b0"]) / theta["b1"] * 
          (exp(theta["b1"] * x) - 1)
      }
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta["c"] * x + exp(theta["b0"]) / theta["b1"] * 
          (exp(theta["b1"] * x) - 1)
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta["a0"]) / theta["a1"] * (1 - exp(-theta["a1"] * x)) + 
          theta["c"] * x + exp(theta["b0"]) / theta["b1"] * 
          (exp(theta["b1"] * x) - 1)
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        (theta["b1"] * x)^theta["b0"]
      }
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta["c"] * x + (theta["b1"] * x)^theta["b0"]
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta["a0"]) / theta["a1"] * (1 - exp(-theta["a1"] * x)) +
          theta["c"] * x + (theta["b1"] * x)^theta["b0"]
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      cumhazfun <- function(theta, x) {
        log(1 + theta["b2"] * exp(theta["b0"]) / theta["b1"] * 
              (exp(theta["b1"] * x) - 1)) * (1 / theta["b2"])
      }
    } else if (algObj$shape == "Makeham") {
      cumhazfun <- function(theta, x) {
        theta["c"] * x + log(1 + theta["b2"] * exp(theta["b0"]) / 
                               theta["b1"] * 
                               (exp(theta["b1"] * x) - 1)) * 
          (1 / theta["b2"])
      }
    } else {
      cumhazfun <- function(theta, x) {
        exp(theta["a0"]) / theta["a1"] * (1 - exp(-theta["a1"] * x)) +
          theta["c"] * x + log(1 + theta["b2"] * 
                                 exp(theta["b0"]) / theta["b1"] * 
                                 (exp(theta["b1"] * x) - 1)) *
          (1 / theta["b2"])
      }
    }
  }
  return(cumhazfun)
}

# c) General survival functions:
.DefineSurvMatrix <- function(algObj) {
  if (algObj$model == "EX") {
    survfun <- function(theta, x) exp(- c(theta) * x)
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        exp(exp(theta[, "b0"]) / theta[, "b1"] * 
              (1 - exp(theta[, "b1"] * x)))
      }
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta[, "c"] * x + exp(theta[, "b0"]) / theta[, "b1"] * 
              (1 - exp(theta[, "b1"] * x)))
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta[, "a0"]) / theta[, "a1"] * (exp(-theta[, "a1"] * x) - 1) - 
              theta[, "c"] * x + exp(theta[, "b0"]) / theta[, "b1"] * 
              (1 - exp(theta[, "b1"] * x)))
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        exp(-(theta[, "b1"] * x)^theta[, "b0"])
      }      
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta[, "c"] * x - (theta[, "b1"] * x)^theta[, "b0"])
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta[, "a0"]) / theta[, "a1"] * (exp(-theta[, "a1"] * x) - 1) -
              theta[, "c"] * x - (theta[, "b1"] * x)^theta[, "b0"])
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        (1 + theta[, "b2"] * exp(theta[, "b0"]) / theta[, "b1"] * 
           (exp(theta[, "b1"] * x) - 1))^(-1 / theta[, "b2"])
      }
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta[, "c"] * x) * (1 + theta[, "b2"] * exp(theta[, "b0"]) / 
                                    theta[, "b1"] * (exp(theta[, "b1"] * x) - 1))^(-1 / theta[, "b2"])
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta[, "a0"]) / theta[, "a1"] * (exp(-theta[, "a1"] * x) - 1) -
              theta[, "c"] * x) * (1 + theta[, "b2"] * 
                                     exp(theta[, "b0"]) / theta[, "b1"] * 
                                     (exp(theta[, "b1"] * x) - 1))^(-1 / theta[, "b2"])
      }
    }
  }
  return(survfun)
}

.DefineSurvNumeric <- function(algObj) {
  if (algObj$model == "EX") {
    survfun <- function(theta, x) exp(- c(theta) * x)
  } else if (algObj$model == "GO") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        exp(exp(theta["b0"]) / theta["b1"] * 
              (1 - exp(theta["b1"] * x)))
      }
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta["c"] * x + exp(theta["b0"]) / theta["b1"] * 
              (1 - exp(theta["b1"] * x)))
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta["a0"]) / theta["a1"] * (exp(-theta["a1"] * x) - 1) - 
              theta["c"] * x + exp(theta["b0"]) / theta["b1"] * 
              (1 - exp(theta["b1"] * x)))
      }
    }
  } else if (algObj$model == "WE") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        exp(-(theta["b1"] * x)^theta["b0"])
      }      
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta["c"] * x - (theta["b1"] * x)^theta["b0"])
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta["a0"]) / theta["a1"] * (exp(-theta["a1"] * x) - 1) -
              theta["c"] * x - (theta["b1"] * x)^theta["b0"])
      }
    }
  } else if (algObj$model == "LO") {
    if (algObj$shape == "simple") {
      survfun <- function(theta, x) {
        (1 + theta["b2"] * exp(theta["b0"]) / theta["b1"] * 
           (exp(theta["b1"] * x) - 1))^(-1 / theta["b2"])
      }
    } else if (algObj$shape == "Makeham") {
      survfun <- function(theta, x) {
        exp(-theta["c"] * x) * (1 + theta["b2"] * exp(theta["b0"]) / 
                                  theta["b1"] * (exp(theta["b1"] * x) - 1))^(-1 / theta["b2"])
      }
    } else {
      survfun <- function(theta, x) {
        exp(exp(theta["a0"]) / theta["a1"] * (exp(-theta["a1"] * x) - 1) -
              theta["c"] * x) * (1 + theta["b2"] * 
                                   exp(theta["b0"]) / theta["b1"] * 
                                   (exp(theta["b1"] * x) - 1))^(-1 / theta["b2"])
      }
    }
  }
  return(survfun)
}

# ------------------------------------------ #
# Likelihood, posterior and prior functions:
# ------------------------------------------ #
# Build initial likelihood object:
.BuildLikeObj <- function(parCovObj, parObj, fullParObj, ageObj,  
                          dataObj, algObj, .CalcMort, .CalcMort.numeric, 
                          .CalcMort.matrix, .CalcSurv, .CalcSurv.numeric, 
                          .CalcSurv.matrix, .CalcCumHaz, .CalcCumHaz.numeric, 
                          .CalcCumHaz.matrix) {
  df <- data.frame(mort = rep(0, dataObj$n), ages = rep(0, dataObj$n))
  likeObj <- list(df = df, pars = c())
  likeObj <- .CalcLike(dataObj, parCovObj, parObj, fullParObj, ageObj, likeObj, 
                       algObj, ind = 1:dataObj$n, .CalcMort, .CalcMort.numeric, 
                       .CalcMort.matrix, .CalcSurv, .CalcSurv.numeric, 
                       .CalcSurv.matrix, .CalcCumHaz, .CalcCumHaz.numeric, 
                       .CalcCumHaz.matrix)
}

# Likelihood function:
.CalcLike <- function(dataObj, parCovObj, parObj, fullParObj, 
                      ageObj, likeObj, algObj, ind, 
                      .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                      .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                      .CalcCumHaz, .CalcCumHaz.numeric, 
                      .CalcCumHaz.matrix) UseMethod(".CalcLike")

# .CalcLike.bastacmr <- function(dataObj, parCovObj, parObj, fullParObj, 
#                                ageObj, likeObj, ind) {
#   # First time created:
#   if (inherits(parObj, "lambda")) {
#     
#     # log-interval censored or survival before min age:
#     likeMbef <- log(exp(- parObj$lambda * ageObj$ages$ageBef[ind]) -
#                       exp(- parObj$lambda * ageObj$ages$ageBef[ind] + 
#                             dataObj$Dx) *
#                       (1 - ageObj$inds$ageBef[ind]))
#     
#     # log-survival truncation before min age:
#     likeHbefTr <- parObj$lambda * ageObj$ages$truBef[ind]
#   } else {
#     likeMbef <- 0
#     likeHbefTr <- 0
#   }
#   # log-mortality after min age:
#   likeMaft <- log(.CalcSurv(parCovObj$theta[ind, ], 
#                             ageObj$ages$ageAft[ind])^exp(parCovObj$gamma[ind]) -
#                     .CalcSurv(parCovObj$theta[ind, ], 
#                               ageObj$ages$ageAft[ind] + 
#                                 dataObj$Dx)^exp(parCovObj$gamma[ind])) * 
#     ageObj$inds$uncens[ind] * (1 - ageObj$inds$ageBef[ind])
#   
#   # log-survival at truncation after min age: 
#   # (CHANGE 2021-01-15, Missed prop. hazards)
#   likeHaftTr <- exp(parCovObj$gamma[ind]) * 
#     .CalcCumHaz(parCovObj$theta[ind, ], ageObj$ages$truAft[ind])
#   
#   # log-likelihood for ages (without recapture):
#   likeAges <- likeMbef + likeMaft
#   
#   # log-likelihood for mortality parameters:
#   likeObj$df$mort[ind] <- likeAges + likeHbefTr + likeHaftTr
#   
#   # Including recapture probability before first obs. and after last obs.:
#   likeRecap <- c((ageObj$alive - dataObj$obsMat)[ind, ] %*% 
#                    log(1 - parObj$pi[fullParObj$pi$idpi]))
#   
#   # log-likelihood for ages:
#   likeObj$df$ages[ind] <- likeAges + likeRecap
#   
#   # Sum of log-like for parameters:
#   likeObj$pars <- sum(likeObj$df$mort)
#   
#   # Assign class to likelihood object:
#   class(likeObj) <- "likecmr"
#   return(likeObj)
# }


.CalcLike.bastacmr <- function(dataObj, parCovObj, parObj, fullParObj, 
                               ageObj, likeObj, algObj, ind, 
                               .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                               .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                               .CalcCumHaz, .CalcCumHaz.numeric, 
                               .CalcCumHaz.matrix) {
  # First time created:
  if (inherits(parObj, "lambda")) {
    
    # log-interval censored or survival before min age:
    likeMbef <- (parObj$lambda * ageObj$ages$ageBef[ind] + dataObj$Dx -
                   parObj$lambda * ageObj$ages$ageBef[ind]) * 
      ageObj$inds$ageBef[ind] -
      parObj$lambda * algObj$minAge * (1 - ageObj$inds$ageBef[ind])
    
    # log-survival truncation before min age:
    likeHbefTr <- parObj$lambda * ageObj$ages$truBef[ind]
  } else {
    likeMbef <- 0
    likeHbefTr <- 0
  }
  
  # log-mortality after min age:
  likeMaft <- log(.CalcSurv(parCovObj$theta[ind, ], 
                            ageObj$ages$ageAft[ind])^exp(parCovObj$gamma[ind]) -
                    .CalcSurv(parCovObj$theta[ind, ], 
                              ageObj$ages$ageAft[ind] + 
                                dataObj$Dx)^exp(parCovObj$gamma[ind])) * 
    ageObj$inds$uncens[ind] * (1 - ageObj$inds$ageBef[ind])
  
  # log-survival at truncation after min age: 
  # (CHANGE 2021-01-15, Missed prop. hazards)
  likeHaftTr <- exp(parCovObj$gamma[ind]) * 
    .CalcCumHaz(parCovObj$theta[ind, ], ageObj$ages$truAft[ind])
  
  # log-likelihood for ages (without recapture):
  likeAges <- likeMbef + likeMaft
  
  # log-likelihood for mortality parameters:
  likeObj$df$mort[ind] <- likeAges + likeHbefTr + likeHaftTr
  
  # Including recapture probability before first obs. and after last obs.:
  likeRecap <- c((ageObj$alive - dataObj$obsMat)[ind, ] %*% 
                   log(1 - parObj$pi[fullParObj$pi$idpi]))
  
  # log-likelihood for ages:
  likeObj$df$ages[ind] <- likeAges + likeRecap
  
  # Sum of log-like for parameters:
  likeObj$pars <- sum(likeObj$df$mort)
  
  # Assign class to likelihood object:
  class(likeObj) <- "likecmr"
  return(likeObj)
}

.CalcLike.bastacensus <- function(dataObj, parCovObj, parObj, fullParObj, 
                                  ageObj, likeObj, algObj, ind, .CalcMort, 
                                  .CalcMort.numeric, .CalcMort.matrix, 
                                  .CalcSurv, .CalcSurv.numeric, 
                                  .CalcSurv.matrix, .CalcCumHaz, 
                                  .CalcCumHaz.numeric, .CalcCumHaz.matrix) {
  # First time created:
  if (inherits(parObj, "lambda")) {
    
    # log-mortality before min age:
    likeMbef <- log(parObj$lambda) * ageObj$inds$ageBef[ind] * 
      ageObj$inds$uncens[ind]
    
    # log-survival before min age:
    likeHbef <- - parObj$lambda * ageObj$ages$ageBef[ind]
    
    # log-survival truncation before min age:
    likeHbefTr <- parObj$lambda * ageObj$ages$truBef[ind]
  } else {
    likeMbef <- 0
    likeHbef <- 0
    likeHbefTr <- 0
  }
  # log-mortality after min age:
  likeMaft <- (parCovObj$gamma[ind] + 
                 log(.CalcMort(parCovObj$theta[ind, ], 
                               ageObj$ages$ageAft[ind]))) * 
    ageObj$inds$uncens[ind] * (1 - ageObj$inds$ageBef[ind])
  
  # log-survival after min age:
  likeHaft <- - exp(parCovObj$gamma[ind]) * 
    .CalcCumHaz(parCovObj$theta[ind, ], ageObj$ages$ageAft[ind])
  
  # log-survival at truncation after min age: 
  # (CHANGE 2021-01-15, Missed prop. hazards)
  likeHaftTr <- exp(parCovObj$gamma[ind]) * 
    .CalcCumHaz(parCovObj$theta[ind, ], ageObj$ages$truAft[ind])
  
  # log-likelihood for ages (without recapture):
  likeAges <- likeMbef + likeHbef + likeMaft + likeHaft
  
  # log-likelihood for mortality parameters:
  likeObj$df$mort[ind] <- likeAges + likeHbefTr + likeHaftTr
  
  # log-likelihood for ages:
  likeObj$df$ages[ind] <- likeAges
  
  # Sum of log-like for parameters:
  likeObj$pars <- sum(likeObj$df$mort)
  
  # Assign class to likelihood object:
  class(likeObj) <- "likecensus"
  return(likeObj)
}


# Create initial posterior object:
.CreatePost <- function(likeObj, fullParObj, parObj, dataObj,
                        ageObj, covObj, priorAgeObj = NULL, .CalcSurv, 
                        .CalcSurv.numeric, .CalcSurv.matrix) {
  # Create posterior list object:
  postObj <- list(priorThe = 0, priorGam = 0, priorLam = 0, postPars = 0, 
                  postAge = rep(0, dataObj$n), priorAge = rep(0, dataObj$n))
  
  # Prior for theta parameters:
  # --------------------------- #
  postObj$priorThe <- sum(.dtnorm(c(parObj$theta), 
                                  c(fullParObj$theta$priorMean), 
                                  c(fullParObj$theta$priorSd), 
                                  lower = c(fullParObj$theta$lower), log = TRUE))
  
  # Prior for gamma parameters:
  # ---------------------------- #
  if (inherits(parObj, "theGam")) {
    postObj$priorGam <- sum(.dtnorm(parObj$gamma, fullParObj$gamma$priorMean,
                                    fullParObj$gamma$priorSd, log = TRUE))
  } 
  
  # Prior for lambda parameter:
  # --------------------------- #
  if (inherits(parObj, "lambda")) {
    postObj$priorLam <- .dtnorm(parObj$lambda, fullParObj$lambda$priorMean,
                                fullParObj$lambda$priorSd, lower = 0, 
                                log = TRUE)
  } 
  
  # Posterior for all parameters:
  # ----------------------------- #
  postObj$postPars <- likeObj$pars + postObj$priorThe + postObj$priorGam + 
    postObj$priorLam
  
  # Posterior for ages:
  # ------------------- #
  if (inherits(dataObj, "bastacmr")) {
    postObj$priorAge <- .CalcPriorAgeDist(ageObj, priorAgeObj, 
                                          ind = 1:dataObj$n, .CalcSurv, 
                                          .CalcSurv.numeric, .CalcSurv.matrix)
  }
  postObj$postAge <- likeObj$df$ages + postObj$priorAge
  
  return(postObj)
}

# Update posterior:
.CalcPost <- function(likeObj, fullParObj, parObj, postObj, dataObj, ageObj,
                      covObj, priorAgeObj = NULL, type = "theta", 
                      ind = NA, .CalcSurv, .CalcSurv.numeric, 
                      .CalcSurv.matrix) {
  if (type == "theta") {
    # Theta parameters:
    # ------------------ #
    postObj$priorThe <- sum(.dtnorm(parObj$theta, fullParObj$theta$priorMean,
                                    fullParObj$theta$priorSd, 
                                    lower = fullParObj$theta$lower, log = TRUE))
  } else if (type == "gamma") {
    # Gamma parameters:
    # ----------------- #
    postObj$priorGam <- sum(.dtnorm(parObj$gamma, fullParObj$gamma$priorMean,
                                    fullParObj$gamma$priorSd, log = TRUE))
  } else if (type == "lambda") {
    # Lambda parameter:
    # ----------------- #
    postObj$priorLam <- .dtnorm(parObj$lambda, fullParObj$lambda$priorMean,
                                fullParObj$lambda$priorSd, lower = 0, 
                                log = TRUE)
  } else {
    # posterior for ages:
    # ------------------- #
    if (inherits(dataObj, "bastacmr")) {
      postObj$priorAge[ind] <- .CalcPriorAgeDist(ageObj, priorAgeObj, 
                                            ind = ind, .CalcSurv, 
                                            .CalcSurv.numeric, 
                                            .CalcSurv.matrix)
    }
  }
  # Posterior for all parameters:
  # ------------------------------ #
  postObj$postPars <- likeObj$pars + postObj$priorThe + postObj$priorGam + 
    postObj$priorLam
  postObj$postAge <- likeObj$df$ages + postObj$priorAge
  return(postObj)
}

# Create object to calculate the prior age distribution:
.SetPriorAgeDist <- function(fullParObj, dataObj, covObj, .CalcSurv, 
                             .CalcSurv.numeric, .CalcSurv.matrix) {
  if (inherits(dataObj, "bastacmr")) {
    dxx <- 0.001
    xx <- seq(0,100,dxx)
    parsPrior <- list()
    parsPrior$theta <- fullParObj$theta$priorMean
    if (inherits(fullParObj, "theGam")) {
      parsPrior$gamma <- fullParObj$gamma$priorMean
    }
    class(parsPrior) <- class(fullParObj)
    priorParsCov <- .CalcCovPars(parsPrior, parCov = NULL, covObj, dataObj, 
                                 type = "both")
    ExVec <- sapply(1:dataObj$n, function(pi) {
      the <- priorParsCov$theta[pi, ]
      gam <- priorParsCov$gamma[pi]
      sum(.CalcSurv(the, xx)^exp(gam) * dxx)
    })
    priorAgeObj <- list(lifeExp = ExVec, pars = parsPrior, 
                        parsCov = priorParsCov)
  } else {
    priorAgeObj <- NA
  }
  return(priorAgeObj)
}

# function to calculate prior age distribution for basta cmr:
.CalcPriorAgeDist <- function(ageObj, priorAgeObj, ind, .CalcSurv,
                              .CalcSurv.numeric, .CalcSurv.matrix) {
  the <- priorAgeObj$parsCov$theta[ind, ]
  gam <- priorAgeObj$parsCov$gamma[ind]
  lifeExp <- priorAgeObj$lifeExp[ind]
  ages <- ageObj$ages$ageAft[ind]
  priorAge <- (.CalcSurv(the, ages)^exp(gam) / lifeExp) *
    (1 - ageObj$inds$ageBef[ind])
  return(priorAge)
}

# Calculate Hastings ratio for Metropolis-Hastings sampling:
.CalcHastRatio <- function(pNow, pNew, parJumps, type = "theta", pp,
                           fullParObj) {
  if (type == 'theta') {
    hRatio <- .dtnorm(pNow$theta[pp], pNew$theta[pp], parJumps$theta[pp], 
                      lower = fullParObj$theta$lower[pp], log = TRUE) -
      .dtnorm(pNew$theta[pp], pNow$theta[pp], parJumps$theta[pp], 
              lower = fullParObj$theta$lower[pp], log = TRUE)
  } else if (type == 'gamma') {
    if (inherits(pNow, "theGam")) {
      hRatio <- .dtnorm(pNow$gamma[pp], pNew$gamma[pp], 
                        parJumps$gamma[pp],
                        lower = fullParObj$gamma$lower[pp], log = TRUE) -
        .dtnorm(pNew$gamma[pp], pNow$gamma[pp], 
                parJumps$gamma[pp],
                lower = fullParObj$gamma$lower[pp], log = TRUE)
    }
  } else {
    if (type == "lambda") {
      hRatio <- .dtnorm(pNow$lambda, pNew$lambda, parJumps$lambda, 
                        lower = fullParObj$lambda$lower, log = TRUE) -
        .dtnorm(pNew$lambda, pNow$lambda, parJumps$lambda, 
                lower = fullParObj$lambda$lower, log = TRUE)
    }
  }
  return(hRatio)
}

# ================================================ #
# ==== D) SAMPLING, JUMP, AND MCMC FUNCTIONS: ====
# ================================================ #
# -------------------- #
# Sampling parameters:
# -------------------- #
# Jitter parameters for different starting values in the 
# MCMC chains:
.JitterPars <- function(parObj, fullParObj, ageObj, dataObj) {
  parObjn <- parObj
  nthe <- fullParObj$theta$len
  parObjn$theta[1:nthe] <- .rtnorm(nthe, c(parObj$theta), 
                                   c(fullParObj$theta$jitter), 
                                   lower = c(fullParObj$theta$lower))
  if (inherits(parObj,  "theGam")) {
    parObjn$gamma <- .rtnorm(fullParObj$gamma$len,
                             fullParObj$gamma$start, 
                             fullParObj$gamma$jitter)
  }
  if (inherits(parObj,  "lambda")) {
    parObjn$lambda <- .rtnorm(fullParObj$lambda$len, 
                              fullParObj$lambda$start, 
                              fullParObj$lambda$jitter, 
                              lower = fullParObj$lambda$lower)
  }
  if (inherits(parObj, "pi")) {
    rho2 <- fullParObj$pi$prior2 + 
      t(t(ageObj$alive - dataObj$Y) %*% rep(1, dataObj$n))
    Rho2 <- tapply(rho2, fullParObj$pi$idpi, sum)
    parObjn$pi <- rbeta(fullParObj$pi$n, fullParObj$pi$Prior1, Rho2)
  }
  return(parObjn)
}

# Sample parameters:
.SamplePars <- function(parObj, parJumps, fullParObj, dataObj,
                        ageObj, type = "theta", pp) {
  parObjn <- parObj
  if (type == 'theta') {
    nthe <- fullParObj$theta$len
    parObjn$theta[pp] <- .rtnorm(1, parObj$theta[pp], parJumps$theta[pp], 
                                 lower = fullParObj$theta$lower[pp])
  } else if (type == 'gamma') {
    if (inherits(parObj, "theGam")) {
      parObjn$gamma[pp] <- .rtnorm(1, parObj$gamma[pp], parJumps$gamma[pp],
                                   lower = fullParObj$gamma$lower[pp])
    }
  } else if (type == "lambda") {
    if (inherits(parObj, "lambda")) {
      parObjn$lambda <- .rtnorm(1, parObj$lambda, parJumps$lambda, 
                                lower = fullParObj$lambda$lower)
    }
  } else {
    if (inherits(parObj, "pi")) {
      rho2 <- fullParObj$pi$prior2 + 
        t(t(ageObj$alive - dataObj$Y) %*% rep(1, dataObj$n))
      Rho2 <- tapply(rho2, fullParObj$pi$idpi, sum)
      parObjn$pi <- rbeta(fullParObj$pi$n, fullParObj$pi$Prior1, Rho2)
    }
  }
  return(parObjn)
}

# -------------- #
# Sampling ages:
# -------------- #
# Sample age object:
.SampleAges <- function(ageObj, dataObj, algObj) UseMethod(".SampleAges")

.SampleAges.agecmr <- function(ageObj, dataObj, algObj) {
  if (dataObj$updA) {
    
    # SAMPLE BIRTH DATES:
    if (dataObj$updB) {
      bi <- ageObj$ages$birth
      
      # Sample births:
      bi[dataObj$idNoB] <- bi[dataObj$idNoB] + 
        sample(-1:1, dataObj$nUpdB, replace = TRUE)
      
      # Find which first obs is not 0:
      idFirst <- dataObj$idNoB[which(dataObj$firstObs[dataObj$idNoB] > 0)]
      # bi[idFirst] <- apply(cbind(bi[idFirst], 
      #                            dataObj$firstObs[idFirst]), 1, min)
      
      if (length(idFirst) > 0) {
        # Find the min between bi and first obs:
        idBiFi <- which(bi[idFirst] > dataObj$firstObs[idFirst] - 1)
        if (length(idBiFi) > 0) {
          bi[idFirst[idBiFi]] <- dataObj$firstObs[idFirst[idBiFi]] - 1
        }
      }
      
      # Find which do not have a first obs:
      idNoFirst <- dataObj$idNoB[which(dataObj$firstObs[dataObj$idNoB] == 0)]
      if (length(idNoFirst) > 0) {
        idBiFi <- which(bi[idNoFirst] > dataObj$ages$death[idNoFirst])
        if (length(idBiFi) > 0) {
          bi[idNoFirst[idBiFi]] <- ageObj$ages$death[idNoFirst[idBiFi]]
        }
        # bi[idNoFirst] <- apply(cbind(bi[idNoFirst], 
        #                              ageObj$ages$death[idNoFirst]), 1, min)
      }
      
      # Find which fall below the minAge:
      if (dataObj$updMinB) {
        idMinB <- dataObj$idMinB[which(bi[dataObj$idMinB] <
                                         dataObj$minBirth[dataObj$idMinB])]
        bi[idMinB] <- dataObj$minBirth[idMinB]
      }
      
      if (dataObj$updMaxB) {
        idMaxB <- dataObj$idMaxB[which(bi[dataObj$idMaxB] >
                                         dataObj$maxBirth[dataObj$idMaxB])]
        bi[idMaxB] <- dataObj$maxBirth[idMaxB]
      }
      
      # Fill in new birth times:
      ageObj$ages$birth <- bi
    }
    
    # SAMPLE DEATH DATES:
    if (dataObj$updD) {
      di <- ageObj$ages$death
      # Sample death year:
      di[dataObj$idNoD] <- di[dataObj$idNoD] +
        sample(-1:1, dataObj$nUpdD, replace = TRUE)
      
      # Find which have a last obs:
      idLast <- dataObj$idNoD[which(dataObj$lastObs[dataObj$idNoD] > 0)]
      if (length(idLast) > 0) {
        # Find which have earlier death than last obs:
        idDiLi <- which(di[idLast] < dataObj$lastObs[idLast])
        if (length(idDiLi) > 0) {
          di[idLast[idDiLi]] <- dataObj$lastObs[idLast[idDiLi]]
        }
      }
      # di[idLast] <- apply(cbind(di[idLast], dataObj$lastObs[idLast]),
      #                     1, max)
      # Find which do no have a last obs:
      idNoLast <- dataObj$idNoD[which(dataObj$lastObs[dataObj$idNoD] == 0)]
      if (length(idNoLast) > 0) {
        # Find which have earlier death than birth:
        idDiLi <- which(di[idNoLast] < ageObj$ages$birth[idNoLast])
        if (length(idDiLi) > 0) {
          di[idNoLast[idDiLi]] <- ageObj$ages$birth[idNoLast[idDiLi]]
        }
      }
      
      # Find which fall below the minAge:
      if (dataObj$updMinD) {
        idMinD <- dataObj$idMinD[which(di[dataObj$idMinD] <
                                         dataObj$minDeath[dataObj$idMinD])]
        di[idMinD] <- dataObj$minDeath[idMinD]
      }
      
      if (dataObj$updMaxD) {
        idMaxD <- dataObj$idMaxD[which(di[dataObj$idMaxD] >
                                         dataObj$maxDeath[dataObj$idMaxD])]
        di[idMaxD] <- dataObj$maxDeath[idMaxD]
      }
      
      # di[idNoLast] <- apply(cbind(di[idNoLast], ageObj$ages$birth[idNoLast]),
      #                       1, max)
      ageObj$ages$death <- di
    }
    ageObj$ages$age <- ageObj$ages$death - ageObj$ages$birth
    ageObj$ages$ageTr <- algObj$start - ageObj$ages$birth
    ageObj$ages$ageTr[which(ageObj$ages$ageTr < 0)] <- 0
    
    # Create alive matrix:
    firstObs <- ageObj$ages$birth + 1
    firstObs[which(firstObs < algObj$start)] <- algObj$start
    # firstObs <- c(apply(cbind(algObj$start, ageObj$ages$birth + 1), 1, max))
    lastObs <- ageObj$ages$death
    # lastObs[which(lastObs > algObj$end)] <- algObj$end
    idlc <- which(lastObs > dataObj$censTime)
    lastObs[idlc] <- dataObj$censTime[idlc]
    # lastObs <- c(apply(cbind(algObj$end, dataObj$censTime, 
    #                          ageObj$ages$death), 1, min))
    ageObj$alive <- .BuildAliveMatrix(firstObs, lastObs, dataObj)
    
    # Assign ages if minAge > 0:
    ageObj <- .AssignMinAge(ageObj, dataObj, algObj)
  }
  return(ageObj)
}

.SampleAges.agecensus <- function(ageObj, dataObj, algObj) {
  if (dataObj$updB) {
    ageObj$ages$birth[dataObj$idNoB] <-     
      round(.rtnorm(dataObj$nUpdB, ageObj$ages$birth[dataObj$idNoB], 
                    0.2, lower = dataObj$bil[dataObj$idNoB], 
                    upper = dataObj$biu[dataObj$idNoB]), 2)
    ageObj$ages$age <- dataObj$lastObs - ageObj$ages$birth
    ageObj$ages$ageTr <- dataObj$firstObs - ageObj$ages$birth
    ageObj$ages$ageTr[ageObj$ages$ageTr < 0] <- 0
    
    # Assign ages if minAge > 0:
    ageObj <- .AssignMinAge(ageObj, dataObj, algObj)
  }
  
  return(ageObj)
}

# Recalculate ages based on minAge:
.AssignMinAge <- function(ageObj, dataObj, algObj) UseMethod(".AssignMinAge")

.AssignMinAge.minAge <- function(ageObj, dataObj, algObj) {
  # Ages after min age:
  ageObj$ages$ageAft <- ageObj$ages$age - algObj$minAge
  ageObj$ages$ageAft[ageObj$ages$ageAft < 0] <- 0
  
  # Ages and indicator before min age:
  ageObj$ages$ageBef <- ageObj$ages$age
  ageObj$inds$ageBef <- rep(0, dataObj$n)
  ageObj$inds$ageBef[which(ageObj$ages$ageBef < algObj$minAge)] <- 1
  ageObj$ages$ageBef[which(ageObj$ages$age >= algObj$minAge)] <- algObj$minAge
  
  # Ages at truncation after min age:
  ageObj$ages$truAft <- ageObj$ages$ageTr - algObj$minAge
  ageObj$ages$truAft[which(ageObj$ages$truAft < 0)] <- 0
  
  # Ages at truncation before min age:
  ageObj$ages$truBef <- ageObj$ages$ageTr
  ageObj$ages$truBef[which(ageObj$ages$ageTr >= algObj$minAge)] <- algObj$minAge
  
  return(ageObj)
}

.AssignMinAge.noMinAge <- function(ageObj, dataObj, algObj) {
  ageObj$ages$ageAft <- ageObj$ages$age
  ageObj$ages$truAft <- ageObj$ages$ageTr
  return(ageObj)
}

# ------------------------------ #
# Dynamic jump update functions:
# ------------------------------ #
# Update jumps:
.UpdateJumps <- function(parJumps, jumpsMat, iter, iterUpd, updTarg) {
  parnames <- names(jumpsMat)
  for (pn in parnames) {
    if (ncol(jumpsMat[[pn]]) > 1) {
      updRate <- apply(jumpsMat[[pn]][iter - ((iterUpd - 1):0), ], 2, sum) / 
        iterUpd  
    } else {
      updRate <- sum(jumpsMat[[pn]][iter - ((iterUpd - 1):0), ]) / 
        iterUpd 
    }
    updRate[updRate == 0] <- 1e-2
    if (is.matrix(parJumps[[pn]])) {
      parJumps[[pn]] <- parJumps[[pn]] * 
        matrix(updRate, nrow(parJumps[[pn]]), ncol(parJumps[[pn]])) / updTarg
    } else {
      parJumps[[pn]] <- parJumps[[pn]] * updRate / updTarg
    }
  }
  return(parJumps)
}

# Function to create parameter jump sd object:
.SetParJumps <- function(parObj, fullParObj) {
  parJump <- parObj
  nthe <- fullParObj$theta$len
  parJump$theta[1:nthe] <- rep(0.1, nthe)
  if (inherits(parObj, "theGam")) {
    parJump$gamma <- rep(0.1, fullParObj$gamma$len)
  }
  if (inherits(parObj, "lambda")) {
    parJump$lambda <- 0.001
  }
  return(parJump)
}

# --------------------------- #
# Functions for MCMC outputs:
# --------------------------- #
# Function to create output and jumps matrices:
.CreateMCMCoutObj <- function(fullParObj, dataObj, algObj, parObj, ageObj,
                              likeObj, postObj, outidObj,
                              type = "mcmc", matrows) {
  McmcOutObj <- list()
  # Create matrices for parameters:
  # matrows <- ifelse(type == "mcmc", algObj$niter, algObj$burnin)
  McmcOutObj$theta <- matrix(0, matrows, fullParObj$theta$len,
                             dimnames = list(NULL, fullParObj$theta$names))
  if (inherits(fullParObj, "theGam")) {
    McmcOutObj$gamma <- matrix(0, matrows, fullParObj$gamma$len, 
                               dimnames = list(NULL, fullParObj$gamma$names))
  }
  if (inherits(fullParObj, "lambda")) {
    McmcOutObj$lambda <- matrix(0, matrows, 1)
  }
  if (inherits(fullParObj, "pi") & type == "mcmc") {
    McmcOutObj$pi <- matrix(0, matrows, fullParObj$pi$len)
  }
  McmcOutObj <- .FillMCMCoutObj(McmcOutObj, parObj, dataObj, ageObj, likeObj,
                                postObj, iter = 1, variables = "params", 
                                type = type)
  # Fill up initial values:
  if (type == "mcmc") {
    # Create matrices for unknown births and deaths:
    if (dataObj$updB) {
      McmcOutObj$B <- matrix(0, outidObj$nKeep, dataObj$nUpdB)
      # McmcOutObj$B[1, ] <- ageObj$ages[dataObj$idNoB, "birth"]
    } else {
      McmcOutObj$B <- NA
    }
    if (inherits(dataObj, "bastacmr")) {
      if (dataObj$updD) {
        McmcOutObj$D <- matrix(0, outidObj$nKeep, dataObj$nUpdD)
        # McmcOutObj$D <- ageObj$ages[dataObj$idNoD, "death"]
      } else {
        McmcOutObj$D <- NA
      }
    } else {
      McmcOutObj$D <- NA
    }
    # Create matrix for likelihood and posterior:
    McmcOutObj$likePost <- matrix(NA, nrow = matrows, ncol = 2,
                                  dimnames = list(NULL, c("Like", "Post")))
    McmcOutObj$likePost[1, ] <- c(likeObj$pars, postObj$postPars)
  }
  return(McmcOutObj)
}

# Function to fill in MCMC output list:
.FillMCMCoutObj <- function(McmcOutObj, parObj, dataObj, ageObj, likeObj,
                            postObj, iter, variables = "params", 
                            type = "mcmc") {
  if (variables == "params") {
    McmcOutObj$theta[iter, ] <- c(parObj$theta)
    if (inherits(parObj, "theGam")) {
      McmcOutObj$gamma[iter, ] <- parObj$gamma    
    }
    if (inherits(parObj, "lambda")) {
      McmcOutObj$lambda[iter, ] <- parObj$lambda
    }
    if (inherits(parObj, "pi") & type == "mcmc") {
      McmcOutObj$pi[iter, ] <- parObj$pi
    }
  } else if (variables == "ages") {
    if (dataObj$updB) {
      McmcOutObj$B[iter, ] <- ageObj$ages[dataObj$idNoB, "birth"]
    }
    if (inherits(dataObj, "bastacmr")) {
      if (dataObj$updD) {
        McmcOutObj$D[iter, ] <- ageObj$ages[dataObj$idNoD, "death"]
      }
    }
  } else {
    McmcOutObj$likePost[iter, ] <- c(likeObj$pars, postObj$postPars)
  }
  return(McmcOutObj)
}

# -------------- #
# MCMC function:
# -------------- #
.RunMCMC <- function(sim, UpdJumps = TRUE, parJumps = NA, ageObj, dataObj,
                     parObj, fullParObj, covObj, priorAgeObj, algObj, outidObj, 
                     .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                     .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                     .CalcCumHaz, .CalcCumHaz.numeric, 
                     .CalcCumHaz.matrix, .JitterPars) {
  # Create initial age object:
  ageNow <- ageObj
  
  # Reset the random seed: 
  runif(sim)
  
  # Create intial parameters:
  parNow <- .JitterPars(parObj, fullParObj, ageObj, dataObj)
  parCovNow <- .CalcCovPars(parObj = parNow, parCov = NULL, 
                            covObj = covObj, dataObj = dataObj, 
                            type = "both")
  likeNow <- .BuildLikeObj(parCovObj = parCovNow, parObj = parNow, 
                           fullParObj = fullParObj, ageObj = ageNow, 
                           dataObj = dataObj, algObj, .CalcMort, 
                           .CalcMort.numeric, .CalcMort.matrix, .CalcSurv,
                           .CalcSurv.numeric, .CalcSurv.matrix, .CalcCumHaz,
                           .CalcCumHaz.numeric, .CalcCumHaz.matrix)
  while (is.na(likeNow$pars) | likeNow$pars == -Inf) {
    parNow <- .JitterPars(parObj, fullParObj, ageObj, dataObj)
    parCovNow <- .CalcCovPars(parObj = parNow, parCov = NULL, 
                              covObj = covObj, dataObj = dataObj, 
                              type = "both")
    likeNow <- .CalcLike(dataObj, parCovNow, parNow, fullParObj, 
                         ageNow, likeObj = NULL, algObj, 
                         ind = 1:dataObj$n,
                         .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                         .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                         .CalcCumHaz, .CalcCumHaz.numeric, .CalcCumHaz.matrix)
  }
  postNow <- .CreatePost(likeObj = likeNow, fullParObj = fullParObj, 
                         parObj = parNow, dataObj = dataObj, ageObj = ageNow, 
                         covObj = covObj, priorAgeObj = priorAgeObj, 
                         .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix)
  if (UpdJumps) {
    # Start jumps for Metropolis algorithm:
    niter <- 7500
    burnin <- niter
    thinning <- 1
    iterUpd <- 50
    updTarg <- 0.25
    updSeq <- seq(1, niter, iterUpd)
    nUpdSeq <- length(updSeq)
    parJumps <- .SetParJumps(parObj, fullParObj)
    jumpsMat <- .CreateMCMCoutObj(fullParObj, dataObj, algObj, 
                                  parObj = parJumps,
                                  ageNow, type = "jumps", matrows = nUpdSeq)
    jumpsUpdMat <- .CreateMCMCoutObj(fullParObj, dataObj, algObj, 
                                     parObj = parJumps,
                                     ageNow, type = "jumps", matrows = niter)
    for (ju in 1:length(jumpsMat)) {
      jumpsUpdMat[[ju]] <- jumpsUpdMat[[ju]] * 0
    }
    
    # Create MCMC output object:
    McmcOutObj <- list()
  } else {
    niter <- algObj$niter
    burnin <- algObj$burnin
    thinning <- algObj$thinning

    # Create MCMC output object:
    McmcOutObj <- .CreateMCMCoutObj(fullParObj, dataObj, algObj, 
                                    parObj = parNow, ageNow, likeNow, 
                                    postNow, outidObj, type = "mcmc", 
                                    matrows = outidObj$nAll)
    
    # Counter for updated ages:
    indAge <- 0
    indPars <- 0
  }
  # Update counter:
  updCount <- parObj
  for (iup in 1:length(updCount)) updCount[[iup]] <- updCount[[iup]] * 0
  
  # Start MCMC:
  for (iter in 1:niter) {
    for (tg in c("theta", 'gamma', "lambda")) {
      if (tg %in% names(fullParObj)) {
        ntg <- fullParObj[[tg]]$len
        for (pp in 1:ntg) {
          parNew <- .SamplePars(parNow, parJumps, fullParObj, dataObj, 
                                ageNow, type = tg, pp = pp)
          parCovNew <- .CalcCovPars(parObj = parNew, parCov = NULL, 
                                    covObj = covObj, dataObj = dataObj, 
                                    type = "both")
          likeNew <- .CalcLike(dataObj, parCovNew, parNew, fullParObj, ageNow, 
                               likeNow, algObj, ind = 1:dataObj$n,
                               .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                               .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                               .CalcCumHaz, .CalcCumHaz.numeric, 
                               .CalcCumHaz.matrix)
          postNew <- .CalcPost(likeNew, fullParObj, parNew, postNow, 
                               dataObj, ageNow, covObj, priorAgeObj, 
                               type = tg, ind = NA, .CalcSurv, 
                               .CalcSurv.numeric, .CalcSurv.matrix)
          hRatio <- .CalcHastRatio(parNow, parNew, parJumps, type = tg,
                                   pp, fullParObj)
          postRatio <- exp(postNew$postPars - postNow$postPars + hRatio)
          if (!is.na(postRatio) & postRatio > runif(1)) {
            parNow <- parNew
            parCovNow <- parCovNew
            likeNow <- likeNew
            postNow <- postNew
            if (iter >= burnin) {
              updCount[[tg]][pp] <- updCount[[tg]][pp] + 1
            }
            if (UpdJumps & iter <= burnin) {
              jumpsUpdMat[[tg]][iter, pp] <- 1
            }
          }
        }
      }
    }

    # Sample recapture probabilities:
    if (inherits(dataObj, "bastacmr")) {
      parNow <- .SamplePars(parNow, parJumps, fullParObj, dataObj, 
                            ageNow, type = "pi", pp = pp)
      if (iter >= burnin) {
        updCount$pi <- updCount$pi + 1
      }
      
      likeNow <- .CalcLike(dataObj, parCovNow, parNow, fullParObj, ageNow, 
                           likeNow, algObj, ind = 1:dataObj$n,
                           .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                           .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                           .CalcCumHaz, .CalcCumHaz.numeric, 
                           .CalcCumHaz.matrix)
      postNow <- .CalcPost(likeNow, fullParObj, parNow, postNow, dataObj, 
                           ageNow, covObj, priorAgeObj, type = "ages", 
                           ind = dataObj$idNoA, .CalcSurv, .CalcSurv.numeric, 
                           .CalcSurv.matrix)
    }
    
    # Sample ages:
    if (dataObj$updA) {
      ageNew <- .SampleAges(ageNow, dataObj, algObj)
      likeNew <- .CalcLike(dataObj, parCovNow, parNow, fullParObj, ageNew, 
                           likeNow, algObj, ind = dataObj$idNoA,
                           .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                           .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                           .CalcCumHaz, .CalcCumHaz.numeric, 
                           .CalcCumHaz.matrix)
      # CHANGE 2021-08-10: Find subset of times of birth and death that 
      #                    did change.
      if (inherits(ageNew, "agecmr")) {
        idNew <- dataObj$idNoA[which(ageNew$ages[dataObj$idNoA, "birth"] != 
                                       ageNow$ages[dataObj$idNoA, "birth"] | 
                                       ageNew$ages[dataObj$idNoA, "death"] != 
                                       ageNow$ages[dataObj$idNoA, "death"])]
      } else {
        idNew <- dataObj$idNoA[which(ageNew$ages[dataObj$idNoA, "birth"] != 
                                       ageNow$ages[dataObj$idNoA, "birth"])]
      }
      
      likeNew <- .CalcLike(dataObj, parCovNow, parNow, fullParObj, ageNew, 
                           likeNow, algObj, ind = idNew,
                           .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                           .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                           .CalcCumHaz, .CalcCumHaz.numeric, .CalcCumHaz.matrix)
      # CHANGE 2021-08-10: Replace ageObj by ageNew:
      postNew <- .CalcPost(likeNew, fullParObj, parNow, postNow, dataObj, 
                           ageNew, covObj, priorAgeObj, type = "age", 
                           ind = idNew, .CalcSurv, .CalcSurv.numeric, 
                           .CalcSurv.matrix)
      postRatio <- exp(postNew$postAge[idNew] - 
                         postNow$postAge[idNew])
      idUpdAges <- idNew[!is.na(postRatio) & 
                           postRatio > runif(length(idNew))]
      
      if (length(idUpdAges) > 0) {
        ageNow$ages[idUpdAges, ] <- ageNew$ages[idUpdAges, ]
        ageNow$inds[idUpdAges, ] <- ageNew$inds[idUpdAges, ]
        if (inherits(dataObj, "bastacmr")) {
          ageNow$alive[idUpdAges, ] <- ageNew$alive[idUpdAges, ]
        }
        # CHANGE 2021-08-10: Missed updating likelihood and posterior:
        likeNow <- .CalcLike(dataObj, parCovNow, parNow, fullParObj, ageNow, 
                             likeNow, algObj, ind = idUpdAges,
                             .CalcMort, .CalcMort.numeric, .CalcMort.matrix, 
                             .CalcSurv, .CalcSurv.numeric, .CalcSurv.matrix,
                             .CalcCumHaz, .CalcCumHaz.numeric, 
                             .CalcCumHaz.matrix)
        postNow <- .CalcPost(likeNow, fullParObj, parNow, postNow, dataObj, 
                             ageNow, covObj, priorAgeObj, type = "age", 
                             ind = idUpdAges, .CalcSurv, .CalcSurv.numeric, 
                             .CalcSurv.matrix)
      }
      
      # Fill up paramters in output object:
      if (!UpdJumps) {
        if (iter %in% outidObj$keep) {
          indAge <- indAge + 1
          McmcOutObj <- .FillMCMCoutObj(McmcOutObj, parNow, dataObj, ageNow,
                                        likeNow, postNow, iter = indAge,
                                        variables = "ages", type = "mcmc")
        }
      }
    }

    # Fill up paramters, likelihood and posterior in output object:
    if (!UpdJumps) {
      if (iter %in% outidObj$all) {
        indPars <- indPars + 1
        McmcOutObj <- .FillMCMCoutObj(McmcOutObj, parNow, dataObj, ageNow,
                                      likeNow, postNow, iter = indPars,
                                      variables = "params",
                                      type = "mcmc")
        McmcOutObj <- .FillMCMCoutObj(McmcOutObj, parNow, dataObj, ageNow,
                                      likeNow, postNow, iter = indPars,
                                      variables = "likepost", type = "mcmc")
        
      }
    }
    
    # Dynamic Metropolis to update jumps sd:
    if (UpdJumps) {
      if (iter %in% updSeq[-1]) {
        idpar <- which(updSeq == iter)
        parJumps <- .UpdateJumps(parJumps, jumpsUpdMat, iter, iterUpd, updTarg)
        # jumpsMat <- .FillMCMCoutObj(jumpsMat, parJumps, dataObj, ageNow, 
        #                             idpar, variables = "params", 
        #                             type = "jumps")
        for (pj in 1:length(jumpsMat)) {
          parName <- names(jumpsMat)[pj]
          jumpsMat[[parName]][idpar, ] <- c(parJumps[[parName]])
        }
        if (iter == max(updSeq)) {
          for (pj in 1:length(jumpsMat)) {
            parName <- names(jumpsMat)[pj]
            len <- fullParObj[[parName]]$len
            njumps <- floor(nUpdSeq / 2):nUpdSeq
            mat <- jumpsMat[[parName]][njumps, ]
            if (is.matrix(mat)) {
              parJumps[[parName]][1:len] <- apply(mat, 2, mean)
            } else {
              parJumps[[parName]][1:len] <- mean(mat)
            }
          }
        }
      }
    }
  }
  McmcOutObj$jumps <- parJumps
  McmcOutObj$update <- updCount
  return(McmcOutObj)
}

# ===================================== #
# ==== F) FUNCTIONS FOR MCMC OUTPUTS: ====
# ===================================== #
# Extract thinned sequences from multiple runs, calculate coefficients,
# DIC and quantiles for mortality, survival, summary statistics and ages:
.ExtractParalOut <- function(bastaOut, outidObj, fullParObj, covObj, covsNames, 
                             nsim, dataObj, algObj, defTheta, .CalcMort, 
                             .CalcMort.numeric, .CalcMort.matrix, 
                             .CalcSurv, .CalcSurv.matrix, 
                             .CalcSurv.numeric) {
  cat("Calculating summary statistics...")
  nthin <- outidObj$nKeep
  parMat <- bastaOut[[1]]$theta[outidObj$idKeep, ]
  # if (inherits(parMat, "numeric")) {
  #   parMat <- matrix(parMat, ncol = 1)
  # }
  parnames <- fullParObj$theta$names
  idTheta <- 1:ncol(bastaOut[[1]]$theta)
  likePost <- bastaOut[[1]]$likePost[outidObj$idKeep, ]
  updVec <- c(bastaOut[[1]]$update$theta)
  updTot <- (algObj$niter - algObj$burnin + 1) * algObj$nsim
  
  # Time of birth estimates:
  birthMat <- matrix(dataObj$bi, dataObj$n, nthin * algObj$nsim)
  
  # Time of death estimates:
  if (algObj$dataType == "census") {
    deathMat <- matrix(dataObj$lastObs, dataObj$n, nthin * algObj$nsim)
  } else {
    deathMat <- matrix(dataObj$di, dataObj$n, nthin * algObj$nsim)
  }
  
  # Parameter matrices:
  if (covsNames$class %in% c("propHaz", "fused")) {
    parMat <- cbind(parMat, bastaOut[[1]]$gamma[outidObj$idKeep, ])
    parnames <- c(parnames, fullParObj$gamma$names)
    updVec <- c(updVec, bastaOut[[1]]$update$gamma)
  }
  if (inherits(fullParObj, "lambda")) {
    parMat <- cbind(parMat, bastaOut[[1]]$lambda[outidObj$idKeep])
    parnames <- c(parnames, "lambda")
    updVec <- c(updVec, bastaOut[[1]]$update$lambda)
  }
  if (inherits(fullParObj, "pi")) {
    parMat <- cbind(parMat, bastaOut[[1]]$pi[outidObj$idKeep, ])
    parnames <- c(parnames, fullParObj$pi$names)
    updVec <- c(updVec, bastaOut[[1]]$update$pi)
  }
  for (sim in 1:nsim) {
    if (dataObj$updB) {
      birthMat[dataObj$idNoB, 1:nthin + (sim - 1) * nthin] <- 
        t(bastaOut[[sim]]$B)
    }
    if (dataObj$updD) {
      deathMat[dataObj$idNoD, 1:nthin + (sim - 1) * nthin] <- 
        t(bastaOut[[sim]]$D)
    }
    if (sim > 1) {
      pmat <- bastaOut[[sim]]$theta[outidObj$idKeep, ]
      tempUpd <- c(bastaOut[[sim]]$update$theta)
      if (inherits(fullParObj, "theGam")) {
        pmat <- cbind(pmat, bastaOut[[sim]]$gamma[outidObj$idKeep, ])
        tempUpd <- c(tempUpd, bastaOut[[sim]]$update$gamma)
      }
      if (inherits(fullParObj, "lambda")) {
        pmat <- cbind(pmat, bastaOut[[sim]]$lambda[outidObj$idKeep])
        tempUpd <- c(tempUpd, bastaOut[[sim]]$update$lambda)
      }
      if (inherits(fullParObj, "pi")) {
        pmat <- cbind(pmat, bastaOut[[sim]]$pi[outidObj$idKeep, ])
        tempUpd <- c(tempUpd, bastaOut[[sim]]$update$pi)
      }
      
      if (inherits(parMat, "matrix")) {
        parMat <- rbind(parMat, pmat)
      } else {
        parMat <- c(parMat, pmat)
        if (sim == nsim) {
          parMat <- matrix(parMat, ncol = 1)
        }
      }
      likePost <- rbind(likePost, bastaOut[[sim]]$likePost[outidObj$idKeep, ])
      updVec <- updVec + tempUpd
    }
  }
  
  # Coefficients for parameters:
  colnames(parMat) <- parnames
  coeffs <- cbind(apply(parMat, 2, mean), apply(parMat, 2, sd), 
                  t(apply(parMat, 2, quantile, c(0.025, 0.975))), 
                  apply(parMat, 2, 
                        function(x) cor(x[-1], x[-length(x)], 
                                        use = "complete.obs")),
                  updVec / updTot)
  colnames(coeffs) <- c("Mean", "StdErr", "Lower95%CI", "Upper95%CI", 
                        "SerAutocorr", "UpdateRate")
  if (nsim > 1) {
    
    idSims <- rep(1:algObj$nsim, each = nthin)
    if (inherits(parMat, "numeric")) {
      parMat <- matrix(parMat, ncol = 1)
    }
    Means <- apply(parMat, 2, function(x) 
      tapply(x, idSims, mean))
    Vars <- apply(parMat, 2, function(x) 
      tapply(x, idSims, var))
    meanall <- apply(Means, 2, mean)
    B <- nthin / (algObj$nsim - 1) * apply(t((t(Means) - meanall)^2), 2, sum)
    W <- 1 / algObj$nsim * apply(Vars, 2, sum)
    Varpl <- (nthin - 1) / nthin * W + 1 / nthin * B
    Rhat <- sqrt(Varpl / W)
    Rhat[which(Varpl == 0)] <- 1
    Rhat[which(Rhat < 1)] <- 1
    conv <- cbind(B, W, Varpl, Rhat)
    rownames(conv) <- colnames(parMat)
    coeffs <- cbind(coeffs, conv[, 'Rhat'])
    colnames(coeffs) <- c(colnames(coeffs)[-ncol(coeffs)], "PotScaleReduc")
    idnconv <- which(conv[, 'Rhat'] > 1.1)
    if (length(idnconv) == 0) {
      # DIC:
      Dave <- mean(- 2 * likePost[, 'Like'])
      pD <- 1/2 * var(-2 * likePost[, 'Like'])
      DIC <- pD + Dave
      Dmode <- Dave - 2 * pD
      k <- fullParObj$theta$len + 
        ifelse("gamma" %in% names(fullParObj), fullParObj$gamma$len, 0) +
        ifelse("lambda" %in% names(fullParObj), fullParObj$lambda$len, 0)
      modSel <- c(Dave, Dmode, pD, k, DIC)
      names(modSel) <- c("D.ave", "D.mode", "pD", "k", "DIC")
      convmessage <- "All parameters converged properly.\n"
    } else {
      modSel <- NA
      convmessage <- "Convergence not reached for some parameters.\n"
      conv <- NULL
    }
  } else {
    modSel <- NA
    convmessage <- "Convergence not calculated due to\ninsuficcient number of simulations.\n"
    conv <- NULL
  }
  cat(" Done.\n")
  cat(convmessage)
  # Kullback-Leibler:
  kulLeib <- .CalcKulbackLeibler(coeffs, covObj, defTheta, fullParObj, 
                                 algObj)
  
  # Agest at death
  ageMat <- deathMat - birthMat
  
  # Age at death or last obs:
  ageLast <- cbind(Mean = apply(ageMat, 1, mean),
                   Median = apply(ageMat, 1, quantile, 0.5),
                   Lower = apply(ageMat, 1, quantile, 0.025),
                   Upper = apply(ageMat, 1, quantile, 0.975))
  
  
  # Age at first detection:
  if (algObj$dataType == "CMR") {
    firstObs <- rep(dataObj$study[1], dataObj$n)
    ageTruncMat <- firstObs - birthMat
    ageTruncMat[which(ageTruncMat < 0)] <- 0
  } else {
    firstObs <- dataObj$firstObs
    ageTruncMat <- firstObs - birthMat
  }
  
  ageFirst <- cbind(Mean = apply(ageTruncMat, 1, mean),
                    Median = apply(ageTruncMat, 1, quantile, 0.5),
                    Lower = apply(ageTruncMat, 1, quantile, 0.025),
                    Upper = apply(ageTruncMat, 1, quantile, 0.975))
  
  # Mortality, survival and density quantiles:
  maxAge <- max(ageLast) * 4
  
  # vector of age intervals for dx values:
  ageInts <- c(0, 1, 5, 10, Inf)
  DxVals <- c(0.001, 0.005, 0.01, 0.05)
  
  # Assign delta x based on maximum longevity:
  Dx <- DxVals[findInterval(maxAge, ageInts)]
  
  # Vector of ages for demographic rate quantiles:
  xv <- seq(0, maxAge, by = Dx)
  
  # Mortality, survival, pdf, and pace-shape quantiles:
  mortQuan <- .CalcDemoFunQuan(parMat, idTheta, xv, covObj, covsNames, defTheta, 
                               funtype = "mort", .CalcMort, 
                               .CalcMort.numeric, .CalcMort.matrix, 
                               .CalcSurv, .CalcSurv.matrix, 
                               .CalcSurv.numeric)
  survQuan <- .CalcDemoFunQuan(parMat, idTheta, xv, covObj, covsNames, defTheta, 
                               funtype = "surv", .CalcMort, 
                               .CalcMort.numeric, .CalcMort.matrix, 
                               .CalcSurv, .CalcSurv.matrix, 
                               .CalcSurv.numeric)
  densQuan <- .CalcDemoFunQuan(parMat, idTheta, xv, covObj, covsNames, defTheta, 
                               funtype = "dens", .CalcMort, 
                               .CalcMort.numeric, .CalcMort.matrix, 
                               .CalcSurv, .CalcSurv.matrix, 
                               .CalcSurv.numeric)
  PSQuan <- .CalcDemoFunQuan(parMat, idTheta, xv, covObj, covsNames, defTheta, 
                             funtype = "PS", .CalcMort, 
                             .CalcMort.numeric, .CalcMort.matrix, 
                             .CalcSurv, .CalcSurv.matrix, 
                             .CalcSurv.numeric)
  # age cuts for basic plots:
  cuts <- list()
  for (nta in names(survQuan)) {
    cuts[[nta]] <- which(survQuan[[nta]][1, ] > 0.02)
  }
  
  # List of outputs:
  outList <- list(params = parMat, coefficients = coeffs, 
                  names = parnames, DIC = modSel, KullbackLeibler = kulLeib, 
                  PS = PSQuan, mort = mortQuan, surv = survQuan, 
                  dens = densQuan, x = xv, cuts = cuts, ageLast = ageLast,
                  ageFirst = ageFirst, convergence = conv, 
                  convmessage = convmessage)
  
  return(outList)
}

# Calculate Kulback-Leibler discrepancies between parameters:
.CalcKulbackLeibler <- function(coeffs, covObj, defTheta, fullParObj, algObj,
                                dataObj) {
  if (!is.null(covObj$cat) & 
      !(length(covObj$cat) == 2 & inherits(covObj, "propHaz"))) {
    if (length(covObj$cat) > 1) {
      if (inherits(covObj, c("fused", "inMort"))) {
        parNames <- defTheta$name
        nPar <- defTheta$length
        lower <- defTheta$lower
        nCat <- length(covObj$cat)
        namesCat <- names(covObj$cat)
      } else {
        parNames <- "gamma"
        nCat <- length(covObj$cat) - 1
        namesCat <- names(covObj$cat)[-1]
        # NOTE 2022-10-14: Fixed issue with categorical covariates on propHaz.
        covCount <- rep(0, length(covObj$genCovs))
        for (icov in 1:length(covObj$genCovs)) {
          idc <- grep(covObj$genCovs[icov], namesCat)
          covCount[icov] <- length(idc)
        }
        nPar <- 1
        lower <- -Inf
      }
      nComb <- (nCat - 1)^2 - ((nCat - 1)^2 - (nCat - 1)) / 2
      covComb1 <- c()
      covComb2 <- c()
      klMat1 <- matrix(0, nPar, nComb, dimnames = list(parNames, NULL))
      klMat2 <- klMat1
      comb <- 0
      for (i in 1:nCat) {
        for (j in 1:nCat) {
          if (i > j) {
            comb <- comb + 1
            covComb1 <- c(covComb1, 
                          sprintf("%s - %s", namesCat[i], namesCat[j]))
            covComb2 <- c(covComb2, 
                          sprintf("%s - %s", namesCat[j], namesCat[i]))
            for (p in 1:nPar) {
              if (inherits(covObj, c("fused", "inMort"))) {
                idP <- sapply(c(i, j), function(ij) 
                  which(rownames(coeffs) == 
                          sprintf("%s.%s", parNames[p], namesCat[ij])))
              } else {
                idP <- sapply(c(i, j), function(ij)
                  grep(namesCat[ij], rownames(coeffs)))
              }
              parRan <- range(sapply(1:2, 
                                     function(pp) .qtnorm(c(0.001, 0.999), 
                                                          coeffs[idP[pp], 1],
                                                          coeffs[idP[pp], 2], 
                                                          lower = lower[p])))
              parVec <- seq(parRan[1], parRan[2], length = 100)
              dp <- parVec[2] - parVec[1]
              parDens <- sapply(1:2, function(pp) 
                .dtnorm(seq(parRan[1], parRan[2], length = 100), 
                        coeffs[idP[pp], 1], coeffs[idP[pp], 2], 
                        lower = lower[p]))
              klMat1[p, comb] <- sum(parDens[, 1] * 
                                       log(parDens[, 1] / parDens[, 2]) * dp)
              klMat2[p, comb] <- sum(parDens[, 2] * 
                                       log(parDens[, 2] / parDens[, 1]) * dp)
            }
          }
        }
      }
      colnames(klMat1) <- covComb1
      colnames(klMat2) <- covComb2
      qKlMat1 <- (1 + (1 - exp(-2 * klMat1)^(1 / 2))) / 2
      qKlMat2 <- (1 + (1 - exp(-2 * klMat2)^(1 / 2))) / 2
      mqKl <- (qKlMat1 + qKlMat2) / 2
      outList <- list(kl1 = klMat1, kl2 = klMat2, qkl1 = qKlMat1, 
                      qkl2 = qKlMat2, mqKl = mqKl)
    } else {
      outList <- "Not calculated"
    }
  } else {
    outList <- "Not calculated"
  }
  return(outList)
}

# Calculate survival and mortality quantiles:
.CalcDemoFunQuan <- function(parMat, idTheta, x, covObj, covsNames, defTheta, 
                             funtype = "mort", .CalcMort, 
                             .CalcMort.numeric, .CalcMort.matrix, 
                             .CalcSurv, .CalcSurv.matrix, 
                             .CalcSurv.numeric) {
  covinf <- list(th = list(), ga = list())
  fullm <- list(th = list(), ga = list())
  nTheta <- length(idTheta)
  if (covsNames$class == "inMort") {
    fullm$th <- parMat[, idTheta]
    covinf$th$num <- ifelse(length(covsNames$cat) == 0, 0, 
                            length(covsNames$cat))
    if (is.na(covsNames$cat)[1]) {
      covinf$th$name <- ""
    } else {
      covinf$th$name <- names(covsNames$cat)
    }
    fullm$ga$cat <- 0
    fullm$ga$con <- 0
    covinf$ga$caname <- ""
    covinf$ga$coname <- ""
    covinf$ga$canum <- 0
    covinf$ga$conum <- 0
  } else if (covsNames$class == "fused") {
    fullm$th <- parMat[, idTheta]
    covinf$th$num <- ifelse(length(covsNames$cat) == 0, 0, 
                            length(covsNames$cat))
    if (is.na(covsNames$cat)[1]) {
      covinf$th$name <- ""
    } else {
      covinf$th$name <- names(covsNames$cat)
    }
    covinf$ga$caname <- ""
    covinf$ga$canum <- 0
    covinf$ga$coname <- names(covsNames$con)
    covinf$ga$conum <- length(covsNames$con)
    fullm$ga$cat <- 0
    concol <- colnames(parMat)[(nTheta + 1):ncol(parMat)]
    idlambda <- grep("lambda", concol)
    if (length(idlambda) > 0) {
      concol <- concol[-idlambda]
    }
    fullm$ga$con <- parMat[, concol]
  } else {
    fullm$th <- parMat[, idTheta]
    covinf$th$num <- 0
    covinf$th$name <- ""
    if (is.na(covsNames$cat)[1]) {
      covinf$ga$caname <- ""
      covinf$ga$canum <- 0
      fullm$ga$cat <- 0
    } else {
      covinf$ga$caname <- names(covsNames$cat)
      covinf$ga$canum <- length(covsNames$cat)
      concol <- colnames(parMat)[(nTheta + 1):ncol(parMat)]
      idlambda <- grep("lambda", concol)
      if (length(idlambda) > 0) {
        concol <- concol[-idlambda]
      }
      idca <- sapply(2:covinf$ga$canum, function(ica) {
        grep(covinf$ga$caname[ica], colnames(parMat))
      })
      fullm$ga$cat <- cbind(0, parMat[, idca])
      colnames(fullm$ga$cat) <- covinf$ga$caname
    }
    if (is.na(covsNames$con)[1]) {
      covinf$ga$coname <- ""
      covinf$ga$conum <- 0
      fullm$ga$con <- 0
    } else {
      covinf$ga$coname <- names(covsNames$con)
      covinf$ga$conum <- length(covsNames$con)
      fullm$ga$con <- parMat[, sprintf("gamma.%s", covinf$ga$coname)]
    }
  }
  idlambda <- grep("lambda", colnames(parMat))
  # if (length(idlambda) == 1) {
  #   etav <- parMat[, idlambda]
  # } else {
  #   etav <- rep(0, nrow(parMat))
  # }
  demoQuan <- list()
  for (nta in covinf$th$name) {
    if (nta == "") {
      ntal <- "nocov"
      catname <- "nocov"
    } else {
      ntal <- nta
      catname <- nta
    }
    #demoQuan[[ntal]] <- list()
    for (nga in covinf$ga$caname) {
      if (nga == "") {
        ngal <- "nocov"
      } else {
        ngal <- nga
        catname <- nga
      }
      if (is.matrix(fullm$th)) {
        thm <- fullm$th[, grep(nta, colnames(fullm$th))]
      } else {
        thm <- matrix(fullm$th, ncol = 1)
      }
      if (!is.matrix(thm)) thm <- matrix(thm, ncol = 1)
      colnames(thm) <- defTheta$name
      if (covinf$ga$canum <= 1) {
        gaca <- fullm$ga$cat
      } else {
        gaca <- fullm$ga$cat[, grep(nga, colnames(fullm$ga$cat))]
      }
      if (covinf$ga$conum == 0) {
        gaco <- fullm$ga$con
      } else {
        if (length(covObj$con) > 1) {
          gaco <- sum(fullm$ga$con * apply(covObj$propHaz[, names(covObj$cont)],
                                           2, mean))
        } else {
          gaco <- fullm$ga$con * mean(covObj$propHaz[, names(covObj$cont)])
        }
      }
      gam <- gaca + gaco
      if (length(gam) == 1) {
        gam <- rep(gam, nrow(thm))
      }
      if (funtype %in% c("mort", "surv", "dens")) {
        demof <- sapply(1:nrow(thm), function(ii) {
          if (funtype == "mort") {
            demf <- .CalcMort(thm[ii, ], x) * exp(gam[ii])
          } else if (funtype == "surv") {
            demf <- .CalcSurv(thm[ii, ], x)^{exp(gam[ii])}
          } else if (funtype == "dens") {
            demf <- .CalcMort(thm[ii, ], x) * exp(gam[ii]) * 
              .CalcSurv(thm[ii, ], x)^{exp(gam[ii])}
          }
          return(demf)
        })
        demofave <- apply(demof, 1, mean, na.rm = TRUE)
        demofci <- apply(demof, 1, quantile, c(0.025, 0.975), na.rm = TRUE)
        demoffin <- rbind(demofave, demofci)
        rownames(demoffin) <- c("Mean", "2.5%", "97.5%")
        demoQuan[[catname]] <- demoffin
      } else {
        Deltax <- x[2] - x[1]
        surv <- sapply(1:nrow(thm), function(ii) {
          sx <- .CalcSurv(thm[ii, ], x)^{exp(gam[ii])}
          return(sx)
        })
        Ex <- apply(surv, 2, .CalcEx, dx = Deltax)
        Hx <- apply(surv, 2, .CalcHx, dx = Deltax)
        Epx <- - log(Hx)
        Gx <- apply(surv, 2, .CalcGx, dx = Deltax)
        PSq <- rbind(c(mean(Ex, na.rm = TRUE), sd(Ex, na.rm = TRUE), 
                       quantile(Ex, c(0.025, 0.975), na.rm = TRUE)),
                     c(mean(Hx, na.rm = TRUE), sd(Hx, na.rm = TRUE), 
                       quantile(Hx, c(0.025, 0.975), na.rm = TRUE)),
                     c(mean(Epx, na.rm = TRUE), sd(Epx, na.rm = TRUE), 
                       quantile(Epx, c(0.025, 0.975), na.rm = TRUE)),
                     c(mean(Gx, na.rm = TRUE), sd(Gx, na.rm = TRUE), 
                       quantile(Gx, c(0.025, 0.975), na.rm = TRUE)))
        dimnames(PSq) <- list(c("LifeExp", "LifeTableEntropy", "LifespanEqual",
                                "Gini"), 
                              c("Mean", "SE", "2.5%", "97.5%"))
        demoQuan[[catname]] <- list(PS = PSq, Ex = Ex, Hx = Hx, Epx = Epx,
                                    Gx = Gx)
      }
    }
  }
  return(demoQuan)
}

# =================================== #
# ==== G) DEMOGRAPHIC FUNCTIONS: ====
# =================================== #
# ------------------------------------------- #
# Functions to calculate pace-shape measures:
# ------------------------------------------- #
# life expectancy:
.CalcEx <- function(Sx, dx) sum(Sx * dx) / Sx[1]

# Keyfitz's entropy:
.CalcHx <- function(Sx, dx) {
  Sx1 <- Sx[Sx > 0]; Sx1 <- Sx1 / Sx1[1]
  -sum(Sx1 * log(Sx1) * dx) / sum(Sx1 * dx)
}

# Gini coefficient:
.CalcGx <- function(Sx, dx) {
  Sx <- Sx / Sx[1]
  Sx <- Sx[Sx > 0]
  return(1 - 1 / sum(Sx * dx) * sum(Sx^2 * dx))
}

# Coefficient of variation:
.CalcCVx <- function(x, Sx, dx) {
  Sx <- Sx / Sx[1]
  idd <- which(Sx > 0)
  Sx <- Sx[idd]
  x <- (x - x[1])[idd]
  dS <- -diff(Sx)
  dS <- dS / sum(dS)
  ex <- sum(Sx * dx)
  return(sqrt(sum((x[-length(x)] + dx/2 - ex)^2 * dS)) / ex)
}

# --------------------- #
# Construct life table:
# --------------------- #
# Function to calculate life table output:
.CalcLifeTable <- function(bastaFinal, covObj, algObj, dataObj) {
  cat("Constructing life table... ")
  
  # Age at death or last obs:
  ageLastQuan <- bastaFinal$ageLast
  ageFirstQuan <- bastaFinal$ageFirst
  
  # Number of levels:
  nq <- ncol(ageFirstQuan)
  
  # level names:
  qnames <- colnames(ageFirstQuan)
  
  # Indicator for censored or death:
  departType <- rep("D", dataObj$n)
  if (algObj$dataType == "census") {
    departType[dataObj$idCens] <- "C"
  }
  
  # Find maximum age:
  maxAge <- max(ageLastQuan)
  
  # vector of age intervals for dx values:
  ageInts <- c(0, 0.5, 1, 3, 10, Inf)
  dxVals <- c(0.05, 0.1, 0.25, 0.5, 1)
  
  # Assign delta x based on maximum longevity:
  dx <- dxVals[findInterval(maxAge, ageInts)]
  
  # Find categorical covariates:
  if (is.null(covObj$cat)) {
    covNames <- c("noCov")
  } else {
    covNames <- names(covObj$cat)
  }
  
  # --------------------- #
  # CALCULATE ESTIMATORS: 
  # --------------------- #
  # Create Life tables:
  LT <- list()
  for (covar in covNames) {
    if (covar == "noCov") {
      idx <- 1:dataObj$n
    } else {
      if (inherits(covObj, c("fused", "inMort"))) {
        covcat <- "inMort"
      } else {
        covcat <- "propHaz"
      }
      if (covcat == "propHaz" & !covar %in% colnames(covObj[[covcat]])) {
        if (length(covNames) > 2) {
          idx <- which(apply(covObj[[covcat]][, covNames[-1]], 1, sum) == 0)
        } else {
          idx <- which(covObj[[covcat]][, covNames[-1]] == 0)
        }
      } else {
        idx <- which(covObj[[covcat]][, covar] == 1)
      }
    }
    # Non-parametric survival output:
    LT[[covar]] <- list()
    
    # Extract mean ages:
    ageLast <- ageLastQuan[idx, 1]
    depType <- departType[idx]
    ageFirst <- ageFirstQuan[idx, 1]
    ageFirst[ageFirst < 0] <- 0
    
    # Non-parametric survival:
    npSx <- .CalcNonParamSurv(ageLast = ageLast, ageFirst = ageFirst, 
                              departType = depType, dx = dx)
    
    
    # Fill up list:
    LT[[covar]]$Mean <- npSx$lt
    LT[[covar]]$ple <- npSx$ple
  }
  cat("done.\n")
  return(LT)
}

# Non-parametric survival estimation:
.CalcNonParamSurv <- function(ageLast, ageFirst = NULL,
                              departType, dx = 1) {
  # Number of records:
  n <- length(ageLast)

  # Avoid minor rounding errors:
  ageLast <- round(ageLast, 8)
  
  # Set age first to 0 if NULL:
  if (is.null(ageFirst)) {
    ageFirst <- rep(0, n)
  } else {
    ageFirst <- round(ageFirst, 8)
  }

  # ------------------------ #
  # Product Limit Estimator:
  # ------------------------ #
  # Find records with same first and last age:
  idsame <- which(ageLast == ageFirst)

  # Increase last age by one day:
  ageLastb <- ageLast
  ageLastb[idsame] <- ageLast[idsame] + 1/365.25

  # Sort ages:
  idsort <- sort.int(ageLastb, index.return = TRUE)$ix

  # Create new age vector (sorted):
  agev <- unique(ageLastb[idsort])

  # Number of ages:
  nage <- length(agev)

  # Cx and delta:
  Cx <- rep(0, nage)
  delx <- rep(0, nage)

  # Fill up Cx and delta:
  for (ii in 1:nage) {
    idNx <- which(ageFirst <= agev[ii] & ageLastb >= agev[ii])
    Cx[ii] <- length(idNx)
    idd <- which(ageLastb == agev[ii] & departType == "D")
    delx[ii] <- length(idd)
  }

  # Calculate product limit estimator:
  ple <- cumprod((1 - 1 / Cx)^delx)

  # Add age 0:
  Ages <- agev
  if (Ages[1] > 0) {
    Ages <- c(0, Ages)
    ple <- c(1, ple)
  }

  # Output:
  pleTab <- data.frame(Ages = Ages, ple = ple, event = c(0, delx))

  # ----------- #
  # Life table:
  # ----------- #
  # Unit age vector for that sex:
  agev <- seq(from = 0, to = ceiling(max(ageLast[which(departType == "D")])),
              by = dx)
  nage <- length(agev)

  # Outputs:
  lx <- Nx <- Dx <- ax <- rep(0, nage)
  lx[1] <- 1
  for (ix in 1:nage) {

    if (ix > 1) {
      idxi <- which(pleTab$Ages > agev[ix])
      if (length(idxi) > 0) {
        lx[ix] <- pleTab$ple[idxi[1] - 1]
      } else {
        lx[ix] <- lx[ix - 1]
      }
    }

    # A) EXPOSURES:
    # Find how many entered the interval (including truncated):
    idNx <- which(ageFirst < agev[ix] + dx & ageLast >= agev[ix])
    nNx <- length(idNx)

    # Extract ages and departType:
    xFirst <- ageFirst[idNx]
    xLast <- ageLast[idNx]
    dType <- departType[idNx]

    # Index for individuals dying within interval:
    idDx <- which(xLast < agev[ix] + dx & dType == "D")
    nDx <- length(idDx)

    # Index of truncated in interval:
    idtr <- which(xFirst >= agev[ix])

    # Index of censored in the interval:
    idce <- which(xLast < agev[ix] + dx & dType == "C")

    # Porportion lived within interval:
    intr <- rep(0, nNx)
    ince <- rep(dx, nNx)
    intr[idtr] <- xFirst[idtr] - agev[ix]
    ince[idce] <- xLast[idce] - agev[ix]
    lived <- (ince - intr) / dx

    # Fill in Nx:
    Nx[ix] <- sum(lived)

    # B) DEATHS:
    # Fill in Dx:
    Dx[ix] <- nDx

    # C) PROPORTION LIVED BY THOSE THAT DIED IN INTERVAL:
    if (Dx[ix] > 1) {
      ylived <- xLast[idDx] - agev[ix]
      ax[ix] <- sum(ylived) / dx / nDx
    } else {
      ax[ix] <- 0
    }
  }

  # Age-specific survival probability:
  px <- rep(0, nage)
  px[-nage] <- lx[-1] / lx[-nage]
  px[which(is.na(px) | px == Inf)] <- 0

  # Age-specific mortality probability:
  qx <- 1 - px

  # Number of individual years lived within the interval:
  Lx <- lx * (1 - ax * qx)
  Lx[is.na(Lx)] <- 0

  # Total number of individual years lived after age x:
  Tx <- rev(cumsum(rev(Lx))) * dx

  # Remaining life expectancy after age x:
  ex <- Tx / lx
  ex[which(is.na(ex))] <- 0

  # Life-table:
  lt <- data.frame(Ages = agev, Nx = Nx, Dx = Dx, lx = lx, px = px,
                   qx = qx, Lx = Lx, Tx = Tx, ex = ex)

  npSurv <- list(ple = pleTab, lt = lt)
  class(npSurv) <- "nonParamSurv"
  return(npSurv)
}

Try the BaSTA package in your browser

Any scripts or data that you put into this service are public.

BaSTA documentation built on June 13, 2025, 9:08 a.m.