R/cdftest_localmean_prop.R

In spsurvey: Spatial Sampling Design and Analysis

################################################################################
# Function: cdftest_localmean_prop (not exported)
# Programmer: Tom Kincaid
# Date: October 22, 2020
# Revised: November 2, 2020 to correctly process the column variable when it
#          includes missing (NA) values
# Revised: May 27 2021 to eliminate use of the finite population correction
#          factor with the local mean variance estimator
#
#' Local Mean Variance/Covarince Estimates of Estimated Population Proportions
#'
#' This function organizes input and output for calculation of the local mean
#' variance/covariance estimator for estimated proportions for categorical data.
#'
#' @param design Object of class \code{survey.design} that specifies a complex
#'   survey design.
#'
#' @param design_names Character vector that provides names of survey design
#'   variables in the design argument.
#'
#' @param warn_ind Logical value that indicates whether warning messages were
#'   generated.
#'
#' @param warn_df Data frame for storing warning messages.
#'
#' @param warn_vec Character vector that contains a subpopulation name, the
#'   first subpopulation level, the second subpopulation level, and an indicator
#'   name.
#'
#' @param subset_local Logical value indicating whether the local mean
#'   neighbor structure is built from the union of the two compared groups
#'   only (\code{TRUE}, the default behavior) or from the full
#'   stratum, with out-of-domain sites zeroed out rather than dropped
#'   (\code{FALSE}). The default value is \code{TRUE}.
#'
#' @return A list containing the following objects:
#'   \itemize{
#'     \item{\code{varest}}{matrix containing the variance/covariance estimates
#'       for the contingency table proportion estimates}
#'     \item{\code{warn_ind}}{logical variable that indicates whether warning
#'       messages were generated}
#'     \item{\code{warn_df}}{data frame for storing warning messages}
#'   }
#'
#' @author Tom Kincaid \email{Kincaid.Tom@@epa.gov}
#'
#' @keywords survey
#'
#' @noRd
################################################################################

cdftest_localmean_prop <- function(design, design_names, warn_ind, warn_df,
                                   warn_vec, subset_local = TRUE) {
  # This supports the CDF comparison test (cont_cdftest()): rowvar
  # identifies which of the groups being compared a site belongs to, and
  # colvar identifies which CDF-threshold bin its response falls in, so
  # rowvar x colvar forms a contingency table of group-by-bin proportions.
  # This function computes the full variance-covariance matrix of that
  # table's cell proportions (needed for the chi-squared-type test
  # statistics in svychisq_localmean()) using the local mean estimator, via
  # cdftestvar_prop(); this parallels cat_localmean_prop() but produces a
  # covariance matrix rather than a vector of independent variances.
  # Per-stratum covariance matrices are combined with the same
  # population-size-share-weighted formula used elsewhere.

  # Assign a value to the function name variable

  fname <- "cdftest_localmean_prop"

  # For variables that exist in the design$variables data frame, assign survey
  # design variables

  # colvar is NA only for a genuinely missing response (item nonresponse),
  # which is always dropped. rowvar is NA for sites outside the two groups
  # being compared; under subset_local = TRUE those are dropped too
  # (original behavior), but under subset_local = FALSE they are kept here
  # so their location can inform the local mean neighbor structure below,
  # and are zeroed out (not dropped) downstream via subpop_ind.
  dframe <- if (subset_local) {
    subset(design$variables, !(is.na(rowvar) | is.na(colvar)))
  } else {
    subset(design$variables, !is.na(colvar))
  }
  for (i in names(design_names)) {
    if (is.null(design_names[[i]])) {
      eval(parse(text = paste0(i, " <- NULL")))
    } else {
      eval(parse(text = paste0(i, " <- dframe[, \"", design_names[[i]], "\"]")))
    }
  }

  # Indicator (1/0) for domain (subpopulation) membership among the rows
  # retained above; only used when subset_local = FALSE, to zero out
  # out-of-domain rows' contribution to the contingency table while still
  # using their location in the neighbor structure

  subpop_ind_full <- as.numeric(!is.na(dframe$rowvar))

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

  cluster_ind <- !is.null(clusterID)

  # Assign values to weight variables

  if (cluster_ind) {
    wgt1 <- dframe$wgt1
    wgt2 <- dframe$wgt2
  } else {
    wgt <- dframe$wgt
  }

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

  stratum_ind <- !is.null(stratumID)

  # For a stratified design, determine whether the subpopulation contains a
  # single stratum

  if (stratum_ind) {
    stratum <- factor(stratumID)
    stratum_levels <- levels(stratum)
    nstrata <- length(stratum_levels)
    if (nstrata == 1) {
      stratum_ind <- FALSE
    }
  }

  # Branch for a stratified sample

  if (stratum_ind) {
    # Calculate values required for weighting strata

    if (cluster_ind) {
      popsize_hat <- tapply(wgt1 * wgt2, stratum, sum)
      sum_popsize_hat <- sum(wgt1 * wgt2)
    } else {
      popsize_hat <- tapply(wgt, stratum, sum)
      sum_popsize_hat <- sum(wgt)
    }

    # Create the varest matrix

    m <- with(dframe, length(levels(rowvar)) * length(levels(colvar)))
    varest <- matrix(0, m, m)
    temp <- paste(
      "interaction(factor(rowvar), factor(colvar))Subpopulation",
      1:2
    )
    colnames_varest <- paste(rep(temp, length(levels(dframe$colvar))),
      rep(levels(dframe$colvar), each = 2),
      sep = "."
    )
    colnames(varest) <- colnames_varest

    # Calculate variance estimates

    for (i in 1:nstrata) {
      # restrict to stratum i by blanking out rowvar (subset_local = TRUE)
      # or colvar (subset_local = FALSE) for every site outside it, so
      # cdftestvar_prop() (which operates on the full design object)
      # effectively sees only this stratum's rows of the contingency table.
      # Blanking colvar instead of rowvar when subset_local = FALSE keeps
      # rowvar's NA-ness meaning only "outside the domain" (rather than
      # also "outside this stratum"), which cdftestvar_prop() relies on to
      # zero out (rather than drop) out-of-domain sites while still using
      # their location in the neighbor structure. The union
      # is.na(rowvar) | is.na(colvar) used to drop rows is unaffected by
      # which of the two carries the "outside this stratum" flag, so this
      # leaves the subset_local = TRUE path unchanged.
      temp <- design_names$stratumID
      tst <- design$variables[, temp] != stratum_levels[i]
      design_temp <- design
      if (subset_local) {
        design_temp$variables$rowvar[tst] <- NA
      } else {
        design_temp$variables$colvar[tst] <- NA
      }
      stratum_i <- stratumID == stratum_levels[i]
      if (cluster_ind) {
        temp <- cdftestvar_prop(
          design_temp, wgt2[stratum_i], xcoord[stratum_i], ycoord[stratum_i],
          stratum_ind, stratum_levels[i], cluster_ind, clusterID[stratum_i],
          wgt1[stratum_i], xcoord1[stratum_i], ycoord1[stratum_i], warn_ind,
          warn_df, warn_vec,
          subset_local = subset_local, subpop_ind = subpop_ind_full[stratum_i]
        )
      } else {
        temp <- cdftestvar_prop(
          design_temp, wgt[stratum_i], xcoord[stratum_i], ycoord[stratum_i],
          stratum_ind, stratum_levels[i], cluster_ind,
          warn_ind = warn_ind,
          warn_df = warn_df, warn_vec = warn_vec,
          subset_local = subset_local, subpop_ind = subpop_ind_full[stratum_i]
        )
      }
      varest_st <- temp$varest
      warn_ind <- temp$warn_ind
      warn_df <- temp$warn_df

      # Add estimates to the varest matrix

      tst <- colnames_varest %in% colnames(varest_st)
      varest[tst, tst] <- varest[tst, tst] +
        ((popsize_hat[i] / sum_popsize_hat)^2) * varest_st

      # End the loop for strata
    }

    # Create the results list

    results <- list(
      varest = varest,
      warn_ind = warn_ind,
      warn_df = warn_df
    )

    # Branch for an unstratified sample
  } else {
    # Calculate the variance/covariance estimates

    if (cluster_ind) {
      results <- cdftestvar_prop(
        design, wgt2, xcoord, ycoord, stratum_ind, NULL, cluster_ind, clusterID,
        wgt1, xcoord1, ycoord1, warn_ind, warn_df, warn_vec,
        subset_local = subset_local, subpop_ind = subpop_ind_full
      )
    } else {
      results <- cdftestvar_prop(
        design, wgt, xcoord, ycoord, stratum_ind, NULL, cluster_ind,
        warn_ind = warn_ind, warn_df = warn_df, warn_vec = warn_vec,
        subset_local = subset_local, subpop_ind = subpop_ind_full
      )
    }
  }

  # Return results

  results
}