NIcalc: Functions to Calculate a Nature Index

Documented in calculateIndex

#' Calculate Index
#'
#' Calculate the Nature Index for a major habitat, or calculate a thematic
#' index.
#'
#' A time series of index values are calculated for each of one or more NIunits.
#' NIunits as well as years included in the time series are specified in
#' \code{x}.
#'
#' Index calculations and associated terminology are explained in the vignette
#' \code{NatureIndexCalculation}.
#'
#' \code{calculateIndex} produces an extensive output for each index value
#' to facilitate further analyses of the results. In addition to a sample of
#' \code{nsim} draws of the index estimate, output includes data describing
#' indicators, ICunits, and BSunits included in the input data set, imputations
#' if present in the data, weights per combination of BSunits and indicators,
#' calculation parameters etc.
#'
#' The vignette \code{objectsInNIcalc} gives a detailed description of output
#' objects generated by \code{calculateIndex} and provides an overview of
#' methods developed for \code{niOutput}, \code{niSeries}, and \code{niValue}
#' classes.
#'
#' @seealso \code{\link{niInput}} and \code{\link{niImputations}} for more
#' detailed description of input data, \code{\link{niOutput}} for more detailed
#' description of class \code{niOutput} objects. \cr The following methods are
#' available within \code{NIcalc} to display and summarize results produced by
#' \code{calculateIndex}: \code{\link{summary.niOutput}},
#' \code{\link{summary.niSeries}}, \code{\link{plot.niSeries}},
#' \code{\link{summaryWeights}}, \code{\link{plot.niValue}}, and
#' \code{\link{plotWeights}}. \cr \code{calculateIndex} calls a series of
#' other functions within NIcalc: \code{\link{calculateBSunitWeights}},
#' \code{\link{calculateNIunitWeights}}, \code{\link{calculateWeights}},
#' \code{\link{sampleObsMat}},\code{\link{scaleObsMat}},
#' \code{\link{indexCalculation}}, and \code{\link{indexCalculationPerBSunit}}.
#'
#' @name calculateIndex
#' @encoding UTF-8
#' @author Bård Pedersen
#'
#' @param x list of class \code{niInput}.
#' @param imputations optional list of class \code{niImputations}.
#' @param awBSunit character, name of variable in \code{x$BSunit} used to calculate
#'   NIunit weights.
#' @param nsim integer - number of bootstrap simulations, default is 1000.
#' @param ... further arguments passed on to \code{\link{calculateBSunitWeights}}
#'   (\code{fids, tgroups, keys, w}), \code{\link{calculateNIunitWeights}}
#'   (\code{awbs, awBSunit}), and \code{\link{scaleObsMat}}
#'   (\code{truncAtRef}).
#'
#' @return \code{calculateIndex} returns a list of class \code{niOutput} where each
#' list element is itself a list of class \code{niSeries} containing the index results
#' for an NIunit. Each element of an \code{niSeries} list again is a list of
#' class \code{niValue}. Objects of class \code{niValue} contain the results from
#' the calculation of a single Nature Index or thematic index value.
#'
#' @examples
#' \dontrun{
#' # Calculation of a thematic index:
#'
#' themeImputes <- imputeData(x = themeData,
#'                                nSim = 1000,
#'                                transConst = 0.01,
#'                                maxit = 20,
#'                                printFlag = TRUE)
#' carnivoreIndex <- calculateIndex(x = themeData,
#'                                  imputations = themeImputes,
#'                                  nsim = 1000,
#'                                  fids = FALSE,
#'                                  tgroups = FALSE,
#'                                  keys = "ignore",
#'                                  w = 0,
#'                                  awbs = TRUE,
#'                                  awBSunit = "Skog")
#' }
#'
#' @export

calculateIndex <- function(x = NULL,
                           imputations = NULL,
                           awBSunit = NULL,
                           nsim = 1000, ...) {


  y <- match.call()

  # Check if all input objects are of the right class.

  if (length(x) == 0) {
    stop("Required argument 'x' missing with no default", call. = TRUE)
  }

  dataSet <- x

  if (!is.niInput(dataSet)) {
    stop(paste("Object '",y[2],"' is not a list of class 'niInput'.",sep=""), call. = TRUE)
  }

  imputationsPresent <- length(imputations) > 0

  if (imputationsPresent & (!is.niImputations(imputations))) {
    stop(paste("Object '",y[3],"' is not a list of class 'niImputations'.",sep=""), call. = TRUE)
  }


  if (length(awBSunit) == 0) {
    stop("Required argument 'awBSunit' missing with no default", call. = TRUE)
  } else if (length(which(names(dataSet$BSunits) == awBSunit)) == 0) {
    stop(paste("No variable with name '", awBSunit,"' in 'x$BSunits'", sep=""), call. = TRUE)
  }
  BSunitWeightVariable <- dataSet$BSunits[[which(names(dataSet$BSunits) == awBSunit)]]

  # Set method parameters for calculation of weights and indices and parameters determining output format.

  setCalculationParameters <- function(fids = TRUE,
                                       tgroups = TRUE,
                                       keys = "specialWeight",
                                       w = 0.5,
                                       awbs = TRUE,
                                       stochastic = "observations",
                                       truncAtRef = TRUE) {

    return(list(fids = fids,
                tgroups = tgroups,
                keys = keys,
                w = w,
                awbs = awbs,
                stochastic = stochastic,
                truncAtRef = truncAtRef))
  }

  param <- setCalculationParameters(...)
  param$imputations <- imputationsPresent
  param$awBSunit <- awBSunit

  # Record number of years, NIunits, ICunits and their names / IDs

  nYears <- length(dataSet$indicatorValues)
  yearNames <- names(dataSet$indicatorValues)

  nNIunits <- dim(dataSet$NIunits)[2]
  NIunitNames <- dimnames(dataSet$NIunits)[[2]]

  ICunitIds <- unique(as.vector(dataSet$ICunits)[!is.na(as.vector(dataSet$ICunits))])
  nICunits <- length(ICunitIds)

  # Determine whether there are observations from all ICunits for each year

  obs <- imps <- allObs <- NULL

  for (i in yearNames) {
    obss <- dataSet$indicatorValues[[i]]$ICunitId[which((dataSet$indicatorValues[[i]]$expectedValue >= 0) &
                                                          (!is.na(dataSet$indicatorValues[[i]]$expectedValue)))]
    obs <- c(obs,list(obss))
    names(obs)[length(obs)] <- i
  }

  if (imputationsPresent) {
    for (i in yearNames) {
      impss <- imputations$identifiers$ICunitId[imputations$identifiers$year == i]
      imps <- c(imps,list(impss))
      names(imps)[length(imps)] <- i
    }
    for (i in yearNames) {
      allObs[i] <- sum(ICunitIds %in% c(obs[[i]],imps[[i]])) == nICunits
    }
  } else {
    for (i in yearNames) {
      allObs[i] <- sum(ICunitIds %in% obs[[i]]) == nICunits
    }
  }

  # Decide whether weights should be calculated for each year

  weightsEachYear <- FALSE
  weightsEachYear <- sum(!allObs) > 0
  param$weights <- ifelse(weightsEachYear,"Recalculated each year",
                          "Identical all years")
  param$nsim <- nsim

  # Check if nsim and number of imputations per missing indicator observation are equal

  if (imputationsPresent) {
    if (dim(imputations$imputations)[2] < nsim) {
      warning(paste("Imputed values per missing indicator observation, ",
                    dim(imputations$imputations)[2],
                    ", are less than the number of simulations ",nsim,".\n",
                    "The ",dim(imputations$imputations)[2],
                    " values will be reused equally frequently in the simulations.",sep=""))
      newDim <- dim(imputations$imputations)[2] * ((nsim%%dim(imputations$imputations)[2] > 0) +
                                                     nsim%/%dim(imputations$imputations)[2])
      antRep <- ((nsim%%dim(imputations$imputations)[2] > 0) + nsim%/%dim(imputations$imputations)[2])
      newMat <- matrix(rep(imputations$imputations,antRep),
                       nrow = dim(imputations$imputations)[1],
                       ncol = newDim)
      for (i in 1:dim(imputations$imputations)[1]) {
        newMat[i,dim(imputations$imputations)[2]:newDim] <- sample(newMat[i,dim(imputations$imputations)[2]:newDim])
      }
      imputations$imputations <- newMat
      remove(newMat,antRep,newDim)
    }
  }

  message(paste("Indices for NIunits '",paste(NIunitNames,collapse = "', '"),
                "'\nand years '",paste(yearNames,collapse = "', '"),
                "' will be calculated.\n",
                "The ",nYears*nNIunits," index distributions will each be based on  ",nsim,
                " simulations.",
                "\nThere are ", nICunits,
                " ICunits with observations in data set '",y[2],"'.",sep=""))

  # Skriv beskjed om vektene
  if (weightsEachYear) {
    message("\nWeights will be recalculated for each year")
  }

  if ("fidelity" %in% names(dataSet$indicators)) {
    fidelityVariable <- which(match(names(dataSet$indicators), table = c("fidelity"), nomatch = 0) != 0)
  } else {
    matchNames <- c("id", "name","keyElement","functionalGroup","functionalGroupId","scalingModel","scalingModelId")
    fidelityVariable <- which(match(names(dataSet$indicators), table = matchNames, nomatch = 0) == 0)
  }

  # Decide how indicator observations and reference values will be treated through the calculations

  refMethod <- ifelse(param$stochastic %in% c("reference","both"),"sample","constant")
  obsMethod <- ifelse(param$stochastic %in% c("observations","both"),"sample","constant")

  if(!weightsEachYear) {

    # If weights are identical for all years, The sets of ICunits and reference values will also be
    # identical for all years.

    # Calculate weights common for all years

    message("\nCalculating weights that are the same for all years ..... ")

    xxx <- calculateBSunitWeights(ICunits = dataSet$ICunits,
                                  indicators = dataSet$indicators$name,
                                  fidelity = dataSet$indicators[[fidelityVariable]],
                                  trophicGroup = dataSet$indicators$functionalGroup,
                                  keyElement = dataSet$indicators$keyElement,
                                  fids = param$fids,
                                  tgroups = param$tgroups,
                                  keys = param$keys,
                                  w = param$w)
    yyy <- calculateNIunitWeights(BSunits = dataSet$BSunits$name,
                                  NIunits = dataSet$NIunits,
                                  awbs = param$awbs,
                                  awBSunit = BSunitWeightVariable)
    NIWeights <- calculateWeights(BSunitWeights = xxx,NIunitWeights = yyy)

    # Sample reference values common for all years

    message("\nSampling reference values ..... ")

    refType <- rep("tradObs",length(dataSet$referenceValues$ICunitId))
    custom <- !is.na(dataSet$referenceValues$customDistributionUUID)
    refType[custom] <- "customObs"
    refmat <- sampleObsMat(ICunitId = dataSet$referenceValues$ICunitId,
                           value = dataSet$referenceValues$expectedValue,
                           distrib = dataSet$referenceValues$distributionFamilyName,
                           mu = dataSet$referenceValues$distParameter1,
                           sig = dataSet$referenceValues$distParameter2,
                           customDistribution = dataSet$referenceValues$customDistribution,
                           imputations = NULL,
                           obsMethod = refMethod,
                           obsType = refType,
                           impYear = NULL,
                           nsim = nsim)

    # If necessary, sample reference values for the calculation of the bbb statistic

    if (refMethod == "sample") {
      refmatbbb <- sampleObsMat(ICunitId = dataSet$referenceValues$ICunitId,
                                value = dataSet$referenceValues$expectedValue,
                                distrib = dataSet$referenceValues$distributionFamilyName,
                                mu = dataSet$referenceValues$distParameter1,
                                sig = dataSet$referenceValues$distParameter2,
                                customDistribution = dataSet$referenceValues$customDistribution,
                                imputations = NULL,
                                obsMethod = "constant",
                                obsType = refType,
                                impYear = NULL,
                                nsim = nsim)
    } else {
      refmatbbb <- refmat
    }

    indices <- bbb <- indicesBSunits <- bbbBSunits <- indicatorIndices <- NULL

    for (i in yearNames) {

      # For each year: sample bootmat

      message(paste("\nSampling and scaling indicator observations from ",i,"..... "))

      obsType <- rep("tradObs",length(c(obs[[i]],imps[[i]])))
      custom <- !is.na(dataSet$indicatorValues[[i]]$customDistributionUUID)
      obsType[custom] <- "customObs"
      obsType[dataSet$indicatorValues[[i]]$ICunitId %in% imps[[i]]] <- "imputations"

      bootmat <- sampleObsMat(ICunitId = dataSet$indicatorValues[[i]]$ICunitId,
                              value = dataSet$indicatorValues[[i]]$expectedValue,
                              distrib = dataSet$indicatorValues[[i]]$distributionFamilyName,
                              mu = dataSet$indicatorValues[[i]]$distParameter1,
                              sig = dataSet$indicatorValues[[i]]$distParameter2,
                              customDistribution = dataSet$indicatorValues[[i]]$customDistribution,
                              imputations = imputations,
                              obsMethod = obsMethod,
                              obsType = obsType,
                              impYear = i,
                              nsim = nsim)

      # For each year: If necessary, sample indicator observations for the calculation of the bbb statistic

      if (obsMethod == "sample") {
        bootmatbbb <- sampleObsMat(ICunitId = dataSet$indicatorValues[[i]]$ICunitId,
                                   value = dataSet$indicatorValues[[i]]$expectedValue,
                                   distrib = dataSet$indicatorValues[[i]]$distributionFamilyName,
                                   mu = dataSet$indicatorValues[[i]]$distParameter1,
                                   sig = dataSet$indicatorValues[[i]]$distParameter2,
                                   customDistribution = dataSet$indicatorValues[[i]]$customDistribution,
                                   imputations = imputations,
                                   obsMethod = "constant",
                                   obsType = obsType,
                                   impYear = i,
                                   nsim = nsim)
      } else {
        bootmatbbb <- bootmat
      }

      # For each year: organize bootmat and refmat so they conform

      refmati <- refmat[as.character(dimnames(bootmat)[[1]]),]
      refmatbbbi <- refmatbbb[as.character(dimnames(bootmat)[[1]]),]

      # For each year: scale indicator observations

      scaleVec <- dataSet$indicatorValues[[i]]$scalingModel
      names(scaleVec) <- dataSet$indicatorValues[[i]]$ICunitId
      scaleVec <- scaleVec[dimnames(bootmat)[[1]]]

      scaledBootmat <- scaleObsMat(bootmat = bootmat,
                                   refmat = refmati,
                                   scalingModel = scaleVec,
                                   truncAtRef = param$truncAtRef)

      # For each year: scale indicator observations for the calculation of the bbb statistic

      scaledBootmatbbb <- scaleObsMat(bootmat = bootmatbbb,
                                      refmat = refmatbbbi,
                                      scalingModel = scaleVec,
                                      truncAtRef = param$truncAtRef)

      # For each year: calculate indices for all NIunits

      indices[[i]] <- indexCalculation(ICunits = dataSet$ICunits,
                                       NIWeights = NIWeights,
                                       scaledBootmat = scaledBootmat,
                                       nsim = nsim)

      bbb[[i]] <- indexCalculation(ICunits = dataSet$ICunits,
                                   NIWeights = NIWeights,
                                   scaledBootmat = scaledBootmatbbb,
                                   nsim = nsim)

      indicesBSunits[[i]] <- indexCalculationPerBSunit(ICunits = dataSet$ICunits,
                                                       BSWeights = xxx,
                                                       scaledBootmat = scaledBootmat,
                                                       nsim = nsim)
      bbbBSunits[[i]] <- indexCalculationPerBSunit(ICunits = dataSet$ICunits,
                                                   BSWeights = xxx,
                                                   scaledBootmat = scaledBootmatbbb,
                                                   nsim = nsim)

    }

    # Assemble results in a list object of class NIoutput

    outputObject <- NULL
    for (j in dimnames(dataSet$NIunits)[[2]]) {
      outputObjectj <- NULL
      for (i in yearNames) {
        indicatorsj <-  dimnames(NIWeights[[j]]) [[2]] [which(colSums(NIWeights[[j]],na.rm = TRUE) > 0)]
        BSunitsj <- dimnames(NIWeights[[j]]) [[1]] [which(rowSums(NIWeights[[j]],na.rm = TRUE) > 0)]
        ICunitsj <- dataSet$ICunits[BSunitsj,indicatorsj]
        ICunitsinj <- unique(ICunitsj[!(is.na(ICunitsj))])
        indicatorDataj <- dataSet$indicators[dataSet$indicators$name %in% indicatorsj,]
        BSunitDataj <- dataSet$BSunits[dataSet$BSunits$name %in% BSunitsj,]
        indexj <- indices[[i]][[j]]
        bbbj <- bbb[[i]][[j]][1]
        indicesBSunitsj <- indicesBSunits[[i]][BSunitsj,]
        bbbBSunitsj <- bbbBSunits[[i]][BSunitsj,1]
        indexWeightsj <- NIWeights[[j]][BSunitsj,indicatorsj]
        BSunitWeightsj <- xxx[BSunitsj,indicatorsj]
        NIunitWeightsj <- yyy[BSunitsj,j]
        imputationsji <- ICunitsinj[ICunitsinj %in% imps[[i]]]
        imputationsj <- imputations$identifiers[(imputations$identifiers$ICunitId %in% imputationsji) &
                                                  (imputations$identifiers$year == i),c(2,1)]
        if (imputationsPresent) {
          nImputationsj <- dim(imputationsj)[1]
        } else {
          nImputationsj <- 0
        }
        metadataj <- c(length(indicatorsj),length(BSunitsj),
                       length(ICunitsinj),nImputationsj)
        names(metadataj) <- c("nIndicators","nBSunits","nICunits","nImputations")

        outputObjectj[[i]] <- niValue(indexArea = j,
                                      call = y,
                                      calculationParameters = param,
                                      metadata = metadataj,
                                      year = as.integer(i),
                                      indicators = indicatorsj,
                                      indicatorData = indicatorDataj,
                                      ICunits = ICunitsinj,
                                      ICunitMatrix = ICunitsj,
                                      imputations = imputationsj,
                                      BSunits = BSunitsj,
                                      BSunitData = BSunitDataj,
                                      BSunitWeights = BSunitWeightsj,
                                      NIunitWeights = NIunitWeightsj,
                                      BSunitIndices = indicesBSunitsj,
                                      BSunitbbb = bbbBSunitsj,
                                      indexWeights = indexWeightsj,
                                      index = indexj,
                                      bbb = bbbj)

      }
      outputObject[[j]] <- niSeries(outputObjectj)
    }
    outputObject <- niOutput(outputObject)

  } else {

    # If weights are not identical for all years, reference values need to be selected for each year.
    # Sample reference values for all ICunits

    message("\nSampling reference values ..... ")

    refType <- rep("tradObs",length(dataSet$referenceValues$ICunitId))
    custom <- !is.na(dataSet$referenceValues$customDistributionUUID)
    refType[custom] <- "customObs"
    refmat <- sampleObsMat(ICunitId = dataSet$referenceValues$ICunitId,
                           value = dataSet$referenceValues$expectedValue,
                           distrib = dataSet$referenceValues$distributionFamilyName,
                           mu = dataSet$referenceValues$distParameter1,
                           sig = dataSet$referenceValues$distParameter2,
                           customDistribution = dataSet$referenceValues$customDistribution,
                           imputations = NULL,
                           obsMethod = refMethod,
                           obsType = refType,
                           impYear = NULL,
                           nsim = nsim)

    # If necessary, sample reference values for the calculation of the bbb statistic

    if (refMethod == "sample") {
      refmatbbb <- sampleObsMat(ICunitId = dataSet$referenceValues$ICunitId,
                                value = dataSet$referenceValues$expectedValue,
                                distrib = dataSet$referenceValues$distributionFamilyName,
                                mu = dataSet$referenceValues$distParameter1,
                                sig = dataSet$referenceValues$distParameter2,
                                customDistribution = dataSet$referenceValues$customDistribution,
                                imputations = NULL,
                                obsMethod = "constant",
                                obsType = refType,
                                impYear = NULL,
                                nsim = nsim)
    } else {
      refmatbbb <- refmat
    }

    # For each year:

    xxx <- yyy <- NIWeights <- NULL

    indices <- bbb <- indicesBSunits <- bbbBSunits <- indicatorIndices <- NULL

    for (i in yearNames) {

      # For each year: select ICunits

      ICunitIdi <- c(obs[[i]],imps[[i]])
      indicatorValuesi <- dataSet$indicatorValues[[i]][dataSet$indicatorValues[[i]]$ICunitId %in% ICunitIdi,]

      # For each year: modify ICunit matrix

      ICunitsi <- dataSet$ICunits
      ICunitsi[!(ICunitsi %in% ICunitIdi)] <- NA

      # For each year: calculate weights

      xxx[[i]] <- calculateBSunitWeights(ICunits = ICunitsi,
                                         indicators = dataSet$indicators$name,
                                         fidelity = dataSet$indicators[[fidelityVariable]],
                                         trophicGroup = dataSet$indicators$functionalGroup,
                                         keyElement = dataSet$indicators$keyElement,
                                         fids = param$fids,
                                         tgroups = param$tgroups,
                                         keys = param$keys,
                                         w = param$w)
      yyy[[i]] <- calculateNIunitWeights(BSunits = dataSet$BSunits$name,
                                         NIunits = dataSet$NIunits,
                                         awbs = param$awbs,
                                         awBSunit = BSunitWeightVariable)
      NIWeights[[i]] <- calculateWeights(BSunitWeights = xxx[[i]],NIunitWeights = yyy[[i]])


      # For each year: sample bootmat

      message(paste("\nSampling and scaling indicator observations from ",i,"..... "))

      obsType <- rep("tradObs",length(ICunitIdi))
      custom <- !is.na(indicatorValuesi$customDistributionUUID)
      obsType[custom] <- "customObs"
      obsType[indicatorValuesi$ICunitId %in% imps[[i]]] <- "imputations"

      bootmati <- sampleObsMat(ICunitId = indicatorValuesi$ICunitId,
                               value = indicatorValuesi$expectedValue,
                               distrib = indicatorValuesi$distributionFamilyName,
                               mu = indicatorValuesi$distParameter1,
                               sig = indicatorValuesi$distParameter2,
                               customDistribution = indicatorValuesi$customDistribution,
                               imputations = imputations,
                               obsMethod = obsMethod,
                               obsType = obsType,
                               impYear = i,
                               nsim = nsim)

      # For each year: If necessary, sample indicator observations for the calculation of the bbb statistic

      if (obsMethod == "sample") {
        bootmatbbbi <- sampleObsMat(ICunitId = indicatorValuesi$ICunitId,
                                    value = indicatorValuesi$expectedValue,
                                    distrib = indicatorValuesi$distributionFamilyName,
                                    mu = indicatorValuesi$distParameter1,
                                    sig = indicatorValuesi$distParameter2,
                                    customDistribution = indicatorValuesi$customDistribution,
                                    imputations = imputations,
                                    obsMethod = "constant",
                                    obsType = obsType,
                                    impYear = i,
                                    nsim = nsim)
      } else {
        bootmatbbbi <- bootmati
      }

      # For each year: select reference values and organize bootmat and refmat so they conform

      refmati <- refmat[as.character(dimnames(bootmati)[[1]]),]
      refmatbbbi <- refmatbbb[as.character(dimnames(bootmati)[[1]]),]

      # For each year: scale indicator observations

      scaleVeci <- indicatorValuesi$scalingModel
      names(scaleVeci) <- indicatorValuesi$ICunitId
      scaleVeci <- scaleVeci[dimnames(bootmati)[[1]]]

      scaledBootmati <- scaleObsMat(bootmat = bootmati,
                                    refmat = refmati,
                                    scalingModel = scaleVeci,
                                    truncAtRef = param$truncAtRef)

      # For each year: scale indicator observations for the calculation of the bbb statistic

      scaledBootmatbbbi <- scaleObsMat(bootmat = bootmatbbbi,
                                       refmat = refmatbbbi,
                                       scalingModel = scaleVeci,
                                       truncAtRef = param$truncAtRef)

      # For each year: calculate indices for all NIunits

      indices[[i]] <- indexCalculation(ICunits = ICunitsi,
                                       NIWeights = NIWeights[[i]],
                                       scaledBootmat = scaledBootmati,
                                       nsim = nsim)

      bbb[[i]] <- indexCalculation(ICunits = ICunitsi,
                                   NIWeights = NIWeights[[i]],
                                   scaledBootmat = scaledBootmatbbbi,
                                   nsim = nsim)

      indicesBSunits[[i]] <- indexCalculationPerBSunit(ICunits = ICunitsi,
                                                       BSWeights = xxx[[i]],
                                                       scaledBootmat = scaledBootmati,
                                                       nsim = nsim)
      bbbBSunits[[i]] <- indexCalculationPerBSunit(ICunits = ICunitsi,
                                                   BSWeights = xxx[[i]],
                                                   scaledBootmat = scaledBootmatbbbi,
                                                   nsim = nsim)

    }

    # Assemble results in a list object of class NIoutput

    outputObject <- NULL
    for (j in dimnames(dataSet$NIunits)[[2]]) {
      outputObjectj <- NULL
      for (i in yearNames) {
        indicatorsj <-  dimnames(NIWeights[[i]][[j]]) [[2]] [which(colSums(NIWeights[[i]][[j]],na.rm = TRUE) > 0)]
        BSunitsj <- dimnames(NIWeights[[i]][[j]]) [[1]] [which(rowSums(NIWeights[[i]][[j]],na.rm = TRUE) > 0)]
        ICunitsj <- dataSet$ICunits[BSunitsj,indicatorsj]
        ICunitsinj <- unique(ICunitsj[!(is.na(ICunitsj))])
        indicatorDataj <- dataSet$indicators[dataSet$indicators$name %in% indicatorsj,]
        BSunitDataj <- dataSet$BSunits[dataSet$BSunits$name %in% BSunitsj,]
        indexj <- indices[[i]][[j]]
        bbbj <- bbb[[i]][[j]][1]
        indicesBSunitsj <- indicesBSunits[[i]][BSunitsj,]
        bbbBSunitsj <- bbbBSunits[[i]][BSunitsj,1]
        indexWeightsj <- NIWeights[[i]][[j]][BSunitsj,indicatorsj]
        BSunitWeightsj <- xxx[[i]][BSunitsj,indicatorsj]
        NIunitWeightsj <- yyy[[i]][BSunitsj,j]
        ICunitsj[is.na(indexWeightsj)] <- NA
        ICunitsinj <- unique(ICunitsj[!(is.na(ICunitsj))])
        imputationsji <- ICunitsinj[ICunitsinj %in% imps[[i]]]
        imputationsj <- imputations$identifiers[(imputations$identifiers$ICunitId %in% imputationsji) &
                                                  (imputations$identifiers$year == i),c(2,1)]
        if (imputationsPresent) {
          nImputationsj <- dim(imputationsj)[1]
        } else {
          nImputationsj <- 0
        }
        metadataj <- c(length(indicatorsj),length(BSunitsj),
                       length(ICunitsinj),nImputationsj)
        names(metadataj) <- c("nIndicators","nBSunits","nICunits","nImputations")

        outputObjectj[[i]] <- niValue(indexArea = j,
                                      call = y,
                                      calculationParameters = param,
                                      metadata = metadataj,
                                      year = as.integer(i),
                                      indicators = indicatorsj,
                                      indicatorData = indicatorDataj,
                                      ICunits = ICunitsinj,
                                      ICunitMatrix = ICunitsj,
                                      imputations = imputationsj,
                                      BSunits = BSunitsj,
                                      BSunitData = BSunitDataj,
                                      BSunitWeights = BSunitWeightsj,
                                      NIunitWeights = NIunitWeightsj,
                                      BSunitIndices = indicesBSunitsj,
                                      BSunitbbb = bbbBSunitsj,
                                      indexWeights = indexWeightsj,
                                      index = indexj,
                                      bbb = bbbj)
      }
      outputObject[[j]] <- niSeries(outputObjectj)
    }
    outputObject <- niOutput(outputObject)
  }

  return(outputObject)

}