spsurvey: Spatial Survey Design and Analysis

Documented in cont.analysis

################################################################################
# Function: cont.analysis
# Programmers: Tony Olsen
#              Tom Kincaid
# Date: July 18, 2002
# Last Revised: August 19, 2014
#'
#' Continuous Variable Data Analysis for Probability Survey Data
#'
#' This function organizes input and output for analysis of continuous data
#' generated by a probability survey.  Input can be either an object of class
#' spsurvey.analysis (see the documentation for function spsurvey.analysis) or
#' through use of the other arguments to this function.
#'
#' @param sites  Data frame consisting of two variables: the first variable
#'   is site IDs, and the second variable is a logical vector indicating which
#'   sites to use in the analysis.  If spsurvey.obj is not provided, then this
#'   argument is required.  The default is NULL.
#'
#' @param subpop Data frame describing sets of populations and subpopulations
#'   for which estimates will be calculated.  The first variable is site IDs.
#'   Each subsequent variable identifies a Type of population, where the
#'   variable name is used to identify Type.  A Type variable identifies each
#'   site with one of the subpopulations of that Type.  If spsurvey.obj is not
#'   provided, then this argument is required.  The default is NULL.
#'
#' @param design Data frame consisting of design variables.  If spsurvey.obj
#'   is not provided, then this argument is required.  The default is NULL.
#'   Variables should be named as follows:
#'     \describe{
#'       \item{siteID}{Vector of site IDs}
#'       \item{wgt}{Vector of weights, which are either the weights for a
#'         single-stage sample or the stage two weights for a two-stage sample}
#'       \item{xcoord}{Vector of x-coordinates for location, which are either
#'         the x-coordinates for a single-stage sample or the stage two
#'         x-coordinates for a two-stage sample}
#'       \item{ycoord}{Vector of y-coordinates for location, which are either
#'         the y-coordinates for a single-stage sample or the stage two
#'         y-coordinates for a two-stage sample}
#'       \item{stratum}{Vector of the stratum codes for each site}
#'       \item{cluster}{Vector of the stage one sampling unit (primary sampling
#'         unit or cluster) codes for each site}
#'       \item{wgt1}{Vector of stage one weights in a two-stage design}
#'       \item{xcoord1}{Vector of the stage one x-coordinates for location in a
#'         two-stage design}
#'       \item{ycoord1}{Vector of the stage one y-coordinates for location in a
#'         two-stage design}
#'       \item{support}{Vector of support values - for a finite resource, the
#'         value one (1) for a for site.  For an extensive resource, the measure
#'         of the sampling unit associated with a site.  Required for
#'         calculation of finite and continuous population 
#'         correction factors.}
#'       \item{swgt}{Vector of size-weights, which is the stage two size-weight
#'         for a two-stage design.}
#'       \item{swgt1}{Vector of stage one size-weights for a two-stage design.}
#'     } 
#'
#' @param data.cont Data frame of continuous response variables.  The first
#'   variable is site IDs.  Subsequent variables are response variables. Missing
#'   data (NA) is allowed.  The default is NULL.
#'
#' @param sigma Measurement error variance.  This variable must be a vector
#'   containing a value for each response variable and must have the names
#'   attribute set to identify the response variable names.  Missing data (NA)
#'   is allowed.  The default is NULL.
#'
#' @param var.sigma Variance of the measurement error variance.  This variable
#'   must be a vector containing a value for each response variable and must
#'   have the names attribute set to identify the response variable names.
#'   Missing data (NA) is allowed.  The default is NULL.
#'
#' @param popsize Known size of the resource, which is used to perform ratio
#'   adjustment to estimators expressed using measurement units for the resource
#'   and to calculate strata proportions for calculating estimates for a
#'   stratified sample.  For a finite resource, this argument is either the
#'   total number of sampling units or the known sum of size-weights.  For an
#'   extensive resource, this argument is the measure of the resource, i.e.,
#'   either known total length for a linear resource or known total area for an
#'   areal resource.  The argument must be in the form of a list containing an
#'   element for each population Type in the subpop data frame, where NULL is a
#'   valid choice for a population Type.  The list must be named using the
#'   column names for the population Types in subpop. If a population Type
#'   doesn't contain subpopulations, then each element of the list is either a
#'   single value for an unstratified sample or a vector containing a value for
#'   each stratum for a stratified sample, where elements of the vector are
#'   named using the stratum codes.  If a population Type contains
#'   subpopulations, then each element of the list is a list containing an
#'   element for each subpopulation, where the list is named using the
#'   subpopulation names.  The element for each subpopulation will be either a
#'   single value for an unstratified sample or a named vector of values for a
#'   stratified sample.  The default is NULL.\cr\cr
#'     Example popsize for a stratified sample:\cr
#'       popsize = list("Pop 1"=c("Stratum 1"=750,\cr
#'                                "Stratum 2"=500,\cr
#'                                "Stratum 3"=250),\cr
#'                      "Pop 2"=list("SubPop 1"=c("Stratum 1"=350,\cr
#'                                                "Stratum 2"=250,\cr
#'                                                "Stratum 3"=150),\cr
#'                                   "SubPop 2"=c("Stratum 1"=250,\cr
#'                                                "Stratum 2"=150,\cr
#'                                                "Stratum 3"=100),\cr
#'                                   "SubPop 3"=c("Stratum 1"=150,\cr
#'                                                "Stratum 2"=150,\cr
#'                                                "Stratum 3"=75)),\cr
#'                      "Pop 3"=NULL)\cr\cr
#'     Example popsize for an unstratified sample:\cr
#'       popsize = list("Pop 1"=1500,\cr
#'                      "Pop 2"=list("SubPop 1"=750,\cr
#'                                   "SubPop 2"=500,\cr
#'                                   "SubPop 3"=375),\cr
#'                      "Pop 3"=NULL)\cr
#'
#' @param popcorrect Logical value that indicates whether finite or continuous
#'   population correction factors should be employed during variance
#'   estimation, where TRUE = use the correction factor and FALSE = do not use
#'   the correction factor.  The default is FALSE.  To employ the correction
#'   factor for a single-stage sample, values must be supplied for argument
#'   pcfsize and for the support variable of the design argument.  To employ the
#'   correction factor for a two-stage sample, values must be supplied for
#'   arguments N.cluster and stage1size, and for the support variable of the
#'   design argument.
#'
#' @param pcfsize Size of the resource, which is required for calculation of
#'   finite and continuous population correction factors for a single-stage
#'   sample. For a stratified sample this argument must be a vector containing a
#'   value for each stratum and must have the names attribute set to identify
#'   the stratum codes.  The default is NULL.
#'
#' @param N.cluster The number of stage one sampling units in the resource,
#'   which is required for calculation of finite and continuous population
#'   correction factors for a two-stage sample.  For a stratified sample this
#'   argument must be a vector containing a value for each stratum and must have
#'   the names attribute set to identify the stratum codes.  The default is
#'   NULL.
#'
#' @param stage1size Size of the stage one sampling units of a two-stage
#'   sample, which is required for calculation of finite and continuous
#'   population correction factors for a two-stage sample and must have the
#'   names attribute set to identify the stage one sampling unit codes.  For a
#'   stratified sample, the names attribute must be set to identify both stratum
#'   codes and stage one sampling unit codes using a convention where the two
#'   codes are separated by the & symbol, e.g., "Stratum 1&Cluster 1". The
#'   default is NULL.
#'
#' @param sizeweight Logical value that indicates whether size-weights should
#'   be used in the analysis, where TRUE = use the size-weights and FALSE = do
#'   not use the size-weights.  The default is FALSE.
#'
#' @param total Logical value that indicates whether the population total
#'   estimate should be included in the output Pct data frame, where TRUE =
#'   include the total estimate and FALSE = do not include the estimate.  The
#'   default is FALSE.
#'
#' @param vartype The choice of variance estimator, where "Local" = local mean
#'   estimator and "SRS" = SRS estimator.  The default is "Local".
#'
#' @param conf Numeric value for the confidence level.  The default is 95.
#'
#' @param pctval  Vector of the set of values at which percentiles are
#'   estimated.  The default set is: {5, 10, 25, 50, 75, 90, 95}.
#'
#' @param spsurvey.obj List of class spsurvey.analysis that was produced by
#'   the function spsurvey.analysis.  Depending on input to that function, some
#'   elements of the list may be NULL.  The default is NULL.
#'
#' @return A list containing either two or four data frames of population
#'   estimates for all combinations of population Types, subpopulations within
#'   Types, and response variables.  The data frames containing deconvoluted CDF
#'   estimates and deconvoluted percentile estimates are only included in the
#'   output list when an input value for measurement error variance is provided
#'   to the function.  CDF estimates are calculated for both proportion and size
#'   of the population.  Standard error estimates and confidence interval
#'   estimates also are calculated.
#'
#'   The four data frames are:
#'   \describe{
#'     \item{\code{CDF}}{data frame containing the CDF estimates}
#'     \item{\code{Pct}}{data frame containing the percentile estimates
#'       plus population mean, standard deviation, and variance estimates}
#'     \item{\code{CDF.D}}{data frame containing the deconvoluted CDF
#'       estimates}
#'     \item{\code{Pct.D}}{data frame containing the deconvoluted
#'       percentile estimates}
#'   }
#'
#'   If an input value for measurement error variance is not provided to the
#'   function, then CDF.D and Pct.D are assigned the value NULL.
#'
#' @section Other Functions Required:
#'   \describe{
#'     \item{\code{\link{dframe.check}}}{check site IDs, the sites data frame,
#'       the subpop data frame, and the data.cat data frame to assure valid
#'       contents and, as necessary, create the sites data frame and the subpop
#'       data frame}
#'     \item{\code{\link{vecprint}}}{takes an input vector and outputs a
#'       character string with line breaks inserted}
#'     \item{\code{\link{uniqueID}}}{creates unique site IDs by appending a
#'       unique number to each occurrence of a site ID}
#'     \item{\code{\link{input.check}}}{check input values for errors,
#'       consistency, and compatibility with analytical functions}
#'     \item{\code{\link{cdf.est}}}{estimate the cumulative distribution
#'       function (CDF) for the proportion (expressed as percent) and the total
#'       of a response variable}
#'     \item{\code{\link{cdf.decon}}}{estimate the deconvoluted CDF for the
#'       proportion and the total of a response variable}
#'     \item{\code{\link{total.est}}}{estimate the population total, mean,
#'       variance, and standard deviation of a response variable}
#'   }
#'
#' @author
#'  Tony Olsen \email{Olsen.Tony@epa.gov}\cr
#'  Tom Kincaid \email{Kincaid.Tom@epa.gov}
#'
#' @keywords survey
#'
#' @examples
#' # Continuous variable example:
#' mysiteID <- paste("Site", 1:100, sep="")
#' mysites <- data.frame(
#'   siteID=mysiteID,
#'   Active=rep(TRUE, 100))
#' mysubpop <- data.frame(
#'   siteID=mysiteID,
#'   All.Sites=rep("All Sites",100),
#'   Resource.Class=rep(c("Good","Poor"), c(55,45)))
#' mydesign <- data.frame(
#'   siteID=mysiteID,
#'   wgt=runif(100, 10, 100),
#'   xcoord=runif(100),
#'   ycoord=runif(100),
#'   stratum=rep(c("Stratum1", "Stratum2"), 50))
#' ContVar <- rnorm(100, 10, 1)
#' mydata.cont <- data.frame(
#'   siteID=mysiteID,
#'   ContVar=ContVar)
#' mypopsize <- list(
#'   All.Sites=c(Stratum1=3500, Stratum2=2000),
#'   Resource.Class=list(Good=c(Stratum1=2500, Stratum2=1500),
#'                       Poor=c(Stratum1=1000, Stratum2=500)))
#' cont.analysis(sites=mysites, subpop=mysubpop, design=mydesign,
#'   data.cont=mydata.cont, popsize=mypopsize)
#'
#' # Include deconvolution estimates:
#' mydata.cont <- data.frame(
#'   siteID=mysiteID,
#'   ContVar=ContVar,
#'   ContVar1=ContVar + rnorm(100, 0, sqrt(0.25)),
#'   ContVar2=ContVar + rnorm(100, 0, sqrt(0.50)))
#' mysigma <- c(ContVar=NA, ContVar1=0.25, ContVar2=0.50)
#' cont.analysis(sites=mysites, subpop=mysubpop[,1:2], design=mydesign,
#'   data.cont=mydata.cont, sigma=mysigma, popsize=mypopsize[1])
#'
#' @export
################################################################################

cont.analysis <- function(sites = NULL, subpop = NULL, design = NULL,
   data.cont = NULL, sigma = NULL, var.sigma = NULL, popsize = NULL,
   popcorrect = FALSE, pcfsize = NULL, N.cluster = NULL, stage1size = NULL,
   sizeweight = FALSE, total = FALSE, vartype = "Local", conf = 95,
   pctval = c(5,10,25,50,75,90,95), spsurvey.obj = NULL) {

# Create a data frame for warning messages

   warn.ind <- FALSE
   warn.df <- NULL
   fname <- "cont.analysis"

# Begin the section when an object of class "spsurvey.analysis" was input to the
# function

   if(inherits(spsurvey.obj, "spsurvey.analysis")) {

# Assign variables from the input list

      sites <- spsurvey.obj$sites
      subpop <- spsurvey.obj$subpop
      design <- spsurvey.obj$design
      data.cont <- spsurvey.obj$data.cont
      sigma <- spsurvey.obj$sigma
      var.sigma <- spsurvey.obj$var.sigma
      stratum.ind <- spsurvey.obj$stratum.ind
      cluster.ind <- spsurvey.obj$cluster.ind
      popsize <- spsurvey.obj$popsize
      pcfactor.ind <- spsurvey.obj$pcfactor.ind
      pcfsize <- spsurvey.obj$pcfsize
      N.cluster <- spsurvey.obj$N.cluster
      stage1size <- spsurvey.obj$stage1size
      swgt.ind <- spsurvey.obj$swgt.ind
      vartype <- spsurvey.obj$vartype
      conf <- spsurvey.obj$conf
      pctval <- spsurvey.obj$pctval

# Assign some required values from the subpop data frame

      ntypes <- dim(subpop)[2]
      typenames <- names(subpop)

# Begin the section when an object of class "spsurvey.analysis" was not input to
# the function

   } else {

# Check that the required data frames have been provided

      if(is.null(design))
         stop("\nThe design data frame must be provided.")
      if(!is.data.frame(design))
         stop("\nThe design argument must be a data frame.")
      if(is.null(data.cont))
         stop("\nThe data.cont data frame must be provided.")

# Check the design data frame for required names

      design.names <- names(design)
      temp <- match(design.names, c("siteID", "wgt", "xcoord", "ycoord",
         "stratum", "cluster", "wgt1", "xcoord1", "ycoord1", "support", "swgt",
         "swgt1"), nomatch=0)
      if(any(temp == 0)) {
         temp.str <- vecprint(design.names[temp == 0])
         stop(paste("\nThe following names used in the design data frame do not match the required names:\n", temp.str))
      }

# Check the sites data frame, the design data frame, the subpop data frame, and
# the data.cont data frame to assure valid contents

      temp <- dframe.check(sites, design, subpop, NULL, data.cont, NULL,
         design.names)
      sites <- temp$sites
      design <- temp$design
      subpop <- temp$subpop
      data.cont <- temp$data.cont

# Assign variables from the design data frame

      siteID <- design$siteID
      wgt <- design$wgt
      xcoord <- design$xcoord
      ycoord <- design$ycoord
      stratum <- design$stratum
      cluster <- design$cluster
      wgt1 <- design$wgt1
      xcoord1 <- design$xcoord1
      ycoord1 <- design$ycoord1
      support <- design$support
      swgt <- design$swgt
      swgt1 <- design$swgt1

# Check site IDs for repeat values and, as necessary, create unique site IDs and
# output a warning message

      temp <- sapply(split(siteID, siteID), length)
      if(any(temp > 1)) {
         warn.ind <- TRUE
         temp.str <- vecprint(names(temp)[temp > 1])
         warn <- paste("The following site ID values occur more than once among the values that were \ninput to the function:\n", temp.str, sep="")
         act <- "Unique site ID values were created.\n"
         warn.df <- rbind(warn.df, data.frame(func=I(fname),
            subpoptype=NA, subpop=NA, indicator=NA, stratum=NA, warning=I(warn),
            action=I(act)))
         siteID <- uniqueID(siteID)
         subpop[,1] <- siteID
         data.cont[,1] <- siteID
      }

# Assign some required values from the subpop data frame

      ntypes <- dim(subpop)[2]
      typenames <- names(subpop)

# Assign a logical value to the indicator variable for a stratified sample

      stratum.ind <- length(unique(stratum)) > 1

# If the sample is stratified, convert stratum to a factor, determine stratum
# levels, and calculate number of strata

      if(stratum.ind) {
         stratum <- factor(stratum)
         stratum.levels <- levels(stratum)
         nstrata <- length(stratum.levels)
      } else {
         stratum.levels <- NULL
         nstrata <- NULL
      }

# Assign a logical value to the indicator variable for a two stage sample

      cluster.ind <- length(unique(cluster)) > 1

# If the sample has two stages, convert cluster to a factor, determine cluster
# levels, and calculate number of clusters

      if(cluster.ind) {
         if(stratum.ind) {
            cluster.in <- cluster
            cluster <- tapply(cluster, stratum, factor)
            cluster.levels <- sapply(cluster, levels, simplify=FALSE)
            ncluster <- sapply(cluster.levels, length)
         } else {
            cluster <- factor(cluster)
            cluster.levels <- levels(cluster)
            ncluster <- length(cluster.levels)
         }
      }

# Ensure that popsize is a list

      if(!is.null(popsize) && !is.list(popsize))
         stop("\nThe popsize argument must be a list")

# If the population correction factor is to be used, ensure that support values
# are provided

      if(popcorrect && is.null(support))
         stop("\nThe logical value that indicates whether finite or continuous population \ncorrection factors should be employed during variance estimation was set to \nTRUE, but support values were not provided in the design data frame.")

# Assign the value of popcorrect to the indicator variable for use of the
# population correction factor

      pcfactor.ind <- popcorrect

# If the sample uses size-weights, ensure that size weights are provided

      if(sizeweight) {
         if(is.null(swgt))
            stop("\nThe logical value that indicates whether size-weights should be employed in the analysis was set to \nTRUE, but size-weights were not provided in the design data frame.")
         if(cluster.ind && is.null(swgt1))
            stop("\nThe sample has two stages and the logical value that indicates whether size- \nweights should be employed in the analysis was set to TRUE, but stage one \nsize-weights were not provided in the design data frame.")
      }

# Assign the value of sizeweight to the indicator variable for use of size
# weights

      swgt.ind <- sizeweight

# Check the vector of measurement error variance values for correct names

      if(!is.null(sigma)) {
         temp.names <- names(data.cont)[-1]
         if(length(sigma) != length(temp.names))
            stop("\nThe vector of measurement error variance values is not the correct length.")
         if(is.null(names(sigma)))
            stop("\nThe vector of measurement error variance values must be named.")
         temp <- match(temp.names, names(sigma), nomatch=0)
         if(any(temp == 0)) {
            temp.str <- vecprint(temp.names[temp == 0])
            stop(paste("\nThe following names for the response variables do not occur among the names for \nthe vector of measurement error variance values:\n", temp.str, sep=""))
         }
         temp <- match(names(sigma), temp.names, nomatch=0)
         if(any(temp == 0)) {
            temp.str <- vecprint(names(sigma)[temp == 0])
            stop(paste("\nThe following names for the vector of measurement error variance values do not \noccur among the names for the response variables:\n", temp.str, sep=""))
         }
         sigma <- sigma[temp]
      }

# Check the vector of values for variance of the measurement error variance for
# correct names

      if(!is.null(var.sigma)) {
         if(length(var.sigma) != length(temp.names))
            stop("\nThe vector of values for variance of the measurement error variance is not the \ncorrect length.")
         if(is.null(names(var.sigma)))
            stop("\nThe vector of values for variance of the measurement error variance must be \nnamed.")
         temp <- match(temp.names, names(var.sigma), nomatch=0)
         if(any(temp == 0)) {
            temp.str <- vecprint(temp.names[temp == 0])
            stop(paste("\nThe following names for the response variables do not occur among the names for \nthe vector of values for variance of the measurement error variance:\n", temp.str, sep=""))
         }
         temp <- match(names(var.sigma), temp.names, nomatch=0)
         if(any(temp == 0)) {
            temp.str <- vecprint(names(var.sigma)[temp == 0])
            stop(paste("\nThe following names for the vector of values for variance of the measurement \nerror variance do not occur among the names for the response variables:\n", temp.str, sep=""))
         }
         var.sigma <- var.sigma[temp]
      }

# Determine the number of response values

      nresp <- dim(design)[1]

# Check for compatibility of input values

      temp <- input.check(nresp, wgt, sigma, var.sigma, xcoord, ycoord,
         stratum.ind, stratum, stratum.levels, nstrata, cluster.ind, cluster,
         cluster.levels, ncluster, wgt1, xcoord1, ycoord1, popsize,
         pcfactor.ind, pcfsize, N.cluster, stage1size, support, swgt.ind, swgt,
         swgt1, vartype, conf, pctval=pctval, subpop=subpop)
      popsize <- temp$popsize
      pcfsize <- temp$pcfsize
      N.cluster <- temp$N.cluster
      stage1size <- temp$stage1size

# If the sample was stratified and had two stages, then reset cluster to its
# input value

      if(stratum.ind && cluster.ind)
         cluster <- cluster.in

# If the sample has two stages, determine whether there are a sufficient number
# of sites in each stage one sampling unit to allow variance calculation

      if(cluster.ind) {
         for(itype in 2:ntypes) {
            temp <- apply(table(cluster, subpop[,itype]) == 1, 2, sum)
            ind <- tapply(cluster, subpop[,itype], function(x)
               length(unique(x)))
            if(any(temp == ind)) {
               temp.str <- vecprint(names(temp)[temp == ind])
               warn.df <<- warn.df
               stop(paste("\nA variance estimate cannot be calculated since all of the stage one sampling \nunits contain a single stage two sampling unit for the following \nsubpopulation(s) of population ", typenames[itype], ":\n", temp.str, "\nEnter the following command to view the warning messages that were generated: \nwarnprnt() \n", sep=""))
            }
            if(any(temp > 0)) {
               temp.str <- vecprint(names(temp)[temp > 0])
               warn <- paste("Since they include one or more stage one sampling units with a single site, \nthe mean of the variance estimates for stage one sampling units with two or \nmore sites will be used as the variance estimate for stage one sampling units \nwith one site for the following subpopulation(s) of population\n", typenames[itype], ":\n", temp.str, sep="")
               act <- "The mean of the variance estimates will be used.\n"
               warn.df <- rbind(warn.df, data.frame(func=I(fname),
                  subpoptype=NA, subpop=NA, indicator=NA, stratum=NA,
                  warning=I(warn), action=I(act)))
            }
         }
      }

# As necessary, assign missing values to the design variables

      if(is.null(xcoord))
         xcoord <- rep(NA, nresp)
      if(is.null(ycoord))
         ycoord <- rep(NA, nresp)
      if(is.null(stratum))
         stratum <- rep(NA, nresp)
      if(is.null(cluster))
         cluster <- rep(NA, nresp)
      if(is.null(wgt1))
         wgt1 <- rep(NA, nresp)
      if(is.null(xcoord1))
         xcoord1 <- rep(NA, nresp)
      if(is.null(ycoord1))
         ycoord1 <- rep(NA, nresp)
      if(is.null(support))
         support <- rep(NA, nresp)
      if(is.null(swgt))
         swgt <- rep(NA, nresp)
      if(is.null(swgt1))
         swgt1 <- rep(NA, nresp)

# Recreate the design data frame

      design <- data.frame(siteID=siteID, wgt=wgt, xcoord=xcoord, ycoord=ycoord,
         stratum=stratum, cluster=cluster, wgt1=wgt1, xcoord1=xcoord1,
         ycoord1=ycoord1, support=support, swgt=swgt, swgt1=swgt1)

# End the section when an object of class "spsurvey.analysis" was not input to
# the function

   }

# Loop through all response variables

   nvar <- dim(data.cont)[2]
   varnames <- names(data.cont)
   nrow1a <- 0
   nrow2a <- 0
   for(ivar in 2:nvar) {

# Find unique values of the response variable

      cdfval <- sort(unique(data.cont[!is.na(data.cont[,ivar]),ivar]))

# Loop through all types of populations

      for(itype in 2:ntypes) {

# Find unique subpopulations of this type of population

         subpopnames <- levels(factor(subpop[,itype]))

# Loop through all subpopulations of this type

         for(isubpop in 1:length(subpopnames)) {

# Select sites in a subpopulation

            subpop.ind <- subpop[,itype] == subpopnames[isubpop]
            subpop.ind[is.na(subpop.ind)] <- FALSE

# Determine whether the subpopulation is empty

            if(all(is.na(data.cont[subpop.ind,ivar]))) {
               warn.ind <- TRUE
               warn <- paste("Subpopulation", subpopnames[isubpop], "of population type", typenames[itype], "for indicator", varnames[ivar], "\ncontains no data.\n")
               act <- "None.\n"
               warn.df <- rbind(warn.df, data.frame(func=I(fname),
                  subpoptype=I(typenames[itype]),
                  subpop=I(subpopnames[isubpop]), indicator=I(varnames[ivar]),
                  stratum=NA,  warning=I(warn), action=I(act)))
               next
            }

# Determine whether the subpopulation contains a single value

            if(sum(!is.na(data.cont[subpop.ind,ivar])) == 1) {
               warn.ind <- TRUE
               warn <- paste("Subpopulation", subpopnames[isubpop], "of population type", typenames[itype], "for indicator", varnames[ivar], "\ncontains a single value.  No analysis was performed.\n")
               act <- "None.\n"
               warn.df <- rbind(warn.df, data.frame(func=I(fname),
                  subpoptype=I(typenames[itype]),
                  subpop=I(subpopnames[isubpop]), indicator=I(varnames[ivar]),
                  stratum=NA,  warning=I(warn), action=I(act)))
               next
            }

# For a stratified sample, remove values from pcfsize, N.cluster, and stage1size
# for strata that do not occur in the subpopulation

            if(stratum.ind) {
               temp.pcfsize <- pcfsize[!is.na(match(names(pcfsize),
                  unique(design[subpop.ind, 5])))]
               temp.N.cluster <- N.cluster[!is.na(match(names(N.cluster),
                  unique(design[subpop.ind, 5])))]
               temp.stage1size <- stage1size[!is.na(match(names(stage1size),
                  unique(design[subpop.ind, 5])))]
            } else {
               temp.pcfsize <- pcfsize
               temp.N.cluster <- N.cluster
               temp.stage1size <- stage1size
           }

# Select values from popsize

            if(is.list(popsize[[itype-1]]))
               temp.popsize <- popsize[[itype-1]][[isubpop]]
            else
               temp.popsize <- popsize[[itype-1]]

# Calculate estimates for the response variable

            warn.vec <- c(typenames[itype], subpopnames[isubpop],
               varnames[ivar])

            temp <- cdf.est(z=data.cont[subpop.ind,ivar],
               wgt=design[subpop.ind,2], x=design[subpop.ind,3],
               y=design[subpop.ind,4], stratum=design[subpop.ind,5],
               cluster=design[subpop.ind,6], wgt1=design[subpop.ind,7],
               x1=design[subpop.ind,8], y1=design[subpop.ind,9],
               popsize=temp.popsize, popcorrect=pcfactor.ind,
               pcfsize=temp.pcfsize, N.cluster=temp.N.cluster,
               stage1size=temp.stage1size, support=design[subpop.ind,10],
               sizeweight=swgt.ind, swgt=design[subpop.ind,11],
               swgt1=design[subpop.ind,12], vartype=vartype, conf=conf,
               cdfval=cdfval, pctval=pctval, check.ind=FALSE, warn.ind=warn.ind,
               warn.df=warn.df, warn.vec=warn.vec)
            temp1.cont <- temp$Results
            warn.ind <- temp$warn.ind
            warn.df <- temp$warn.df

            if(!is.null(sigma) && !is.na(sigma[ivar-1])) {
               temp.vs <- NULL
               if(!is.null(var.sigma) && !is.na(var.sigma[ivar-1]))
                  temp.vs <- var.sigma[ivar-1]
               temp <- cdf.decon(z=data.cont[subpop.ind,ivar],
                  wgt=design[subpop.ind,2], sigma=sigma[ivar-1],
                  var.sigma=temp.vs, x=design[subpop.ind,3],
                  y=design[subpop.ind,4], stratum=design[subpop.ind,5],
                  cluster=design[subpop.ind,6], wgt1=design[subpop.ind,7],
                  x1=design[subpop.ind,8],y1=design[subpop.ind,9],
                  popsize=temp.popsize, popcorrect=pcfactor.ind,
                  pcfsize=temp.pcfsize, N.cluster=temp.N.cluster,
                  stage1size=temp.stage1size, support=design[subpop.ind,10],
                  sizeweight=swgt.ind, swgt=design[subpop.ind,11],
                  swgt1=design[subpop.ind,12], vartype=vartype, conf=conf,
                  cdfval=cdfval, pctval=pctval, check.ind=FALSE,
                  warn.ind=warn.ind, warn.df=warn.df, warn.vec=warn.vec)
               temp2.cont <- temp$Results
               warn.ind <- temp$warn.ind
               warn.df <- temp$warn.df
            }

            temp <- total.est(z=data.cont[subpop.ind,ivar],
               wgt=design[subpop.ind,2], x=design[subpop.ind,3],
               y=design[subpop.ind,4], stratum=design[subpop.ind,5],
               cluster=design[subpop.ind,6], wgt1=design[subpop.ind,7],
               x1=design[subpop.ind,8], y1=design[subpop.ind,9],
               popsize=temp.popsize, popcorrect=pcfactor.ind,
               pcfsize=temp.pcfsize, N.cluster=temp.N.cluster,
               stage1size=temp.stage1size, support=design[subpop.ind,10],
               sizeweight=swgt.ind, swgt=design[subpop.ind,11],
               swgt1=design[subpop.ind,12], vartype=vartype, conf=conf,
               check.ind=FALSE, warn.ind=warn.ind, warn.df=warn.df,
               warn.vec=warn.vec)
            if(total) {
               temp3.cont <- temp$Results
            } else {
               temp3.cont <- temp$Results[-1,]
            }
            warn.ind <- temp$warn.ind
            warn.df <- temp$warn.df

# Assign estimates for the response variable to data frames

            if(nrow1a == 0) {
               nn <- dim(temp1.cont$CDF)[1]
               cdfsum <- data.frame(Type=rep(typenames[itype],nn),
                  Subpopulation=rep(subpopnames[isubpop],nn),
                  Indicator=rep(varnames[ivar],nn), temp1.cont$CDF)
               nrow1a <- nn

               nn <- dim(temp1.cont$Pct)[1] + dim(temp3.cont)[1]
               pctsum <- data.frame(Type=rep(typenames[itype],nn),
                  Subpopulation=rep(subpopnames[isubpop],nn),
                  Indicator=rep(varnames[ivar],nn),
                  rbind(temp1.cont$Pct[,1:6], temp3.cont))
               nrow1b <- nn
            } else {
               nn <- dim(temp1.cont$CDF)[1]
               cdfsum <- rbind(cdfsum, data.frame(Type=rep(typenames[itype],nn),
                  Subpopulation=rep(subpopnames[isubpop],nn),
                  Indicator=rep(varnames[ivar],nn), temp1.cont$CDF,
                  row.names=(nrow1a+1):(nrow1a+nn)))
               nrow1a <- nrow1a + nn

               nn <- dim(temp1.cont$Pct)[1] + dim(temp3.cont)[1]
               pctsum <- rbind(pctsum, data.frame(Type=rep(typenames[itype],
                  nn), Subpopulation=rep(subpopnames[isubpop],nn),
                  Indicator=rep(varnames[ivar],nn),
                  rbind(temp1.cont$Pct[,1:6], temp3.cont),
                  row.names=(nrow1b+1):(nrow1b+nn)))
               nrow1b <- nrow1b + nn
            }

            if(!is.null(sigma) && !is.na(sigma[ivar-1])) {
               if(nrow2a == 0) {
                  nn <- dim(temp2.cont$CDF)[1]
                  cdfsum.D <- data.frame(Type=rep(typenames[itype],nn),
                     Subpopulation=rep(subpopnames[isubpop],nn),
                     Indicator=rep(varnames[ivar],nn), temp2.cont$CDF)
                  nrow2a <- nn

                  nn <- dim(temp2.cont$Pct)[1]
                  pctsum.D <- data.frame(Type=rep(typenames[itype],nn),
                     Subpopulation=rep(subpopnames[isubpop],nn),
                     Indicator=rep(varnames[ivar],nn), temp2.cont$Pct[,1:6])
                  nrow2b <- nn
               } else {
                  nn <- dim(temp2.cont$CDF)[1]
                  cdfsum.D <- rbind(cdfsum.D, data.frame(Type=rep(typenames
                     [itype], nn), Subpopulation=rep(subpopnames[isubpop],nn),
                     Indicator=rep(varnames[ivar],nn), temp2.cont$CDF,
                     row.names=(nrow2a+1):(nrow2a+nn)))
                  nrow2a <- nrow2a + nn

                  nn <- dim(temp2.cont$Pct)[1]
                  pctsum.D <- rbind(pctsum.D, data.frame(Type=rep(typenames
                     [itype],nn), Subpopulation=rep(subpopnames[isubpop],nn),
                     Indicator=rep(varnames[ivar],nn), temp2.cont$Pct[,1:6],
                     row.names=(nrow2b+1):(nrow2b+nn)))
                  nrow2b <- nrow2b + nn
               }
            }

# End of the loop for subpopulations

         }

# End of the loop for type of population

      }

# End of the loop for response variables

   }

# As necessary, output a message indicating that warning messages were generated
# during execution of the program

   if(warn.ind) {
      warn.df <<- warn.df
      if(nrow(warn.df) == 1)
         cat("During execution of the program, a warning message was generated.  The warning \nmessage is stored in a data frame named 'warn.df'.  Enter the following command \nto view the warning message: warnprnt()\n")
      else
         cat(paste("During execution of the program,", nrow(warn.df), "warning messages were generated.  The warning \nmessages are stored in a data frame named 'warn.df'.  Enter the following \ncommand to view the warning messages: warnprnt() \nTo view a subset of the warning messages (say, messages number 1, 3, and 5), \nenter the following command: warnprnt(m=c(1,3,5))\n"))
   }

# Assign consecutive numbers to the row names of the output data frames

   if(!is.null(sigma)) {
      row.names(cdfsum) <- 1:nrow(cdfsum)
      row.names(pctsum) <- 1:nrow(pctsum)
      row.names(cdfsum.D) <- 1:nrow(cdfsum.D)
      row.names(pctsum.D) <- 1:nrow(pctsum.D)
   } else {
      row.names(cdfsum) <- 1:nrow(cdfsum)
      row.names(pctsum) <- 1:nrow(pctsum)
   }

# Return the data frames as a list

   if(!is.null(sigma))
      list(CDF=cdfsum, Pct=pctsum, CDF.D=cdfsum.D, Pct.D=pctsum.D)
   else
      list(CDF=cdfsum, Pct=pctsum, CDF.D=NULL, Pct.D=NULL)
}