R/gap.R

Defines functions addALPrefix getALNames gmatchAttr pt2 dofCompute parseVars prepGapDf print.gapList roundLikePrintCoefmat print.gap calLinkingError calLinkingErrorPISA gapHelper gap

Documented in gap print.gap print.gapList

#' @title Gap Analysis
#'
#' @description Compares the average levels of a variable between two groups
#'              that potentially share members.
#'
#' @param variable a character indicating the variable to be compared,
#'                 potentially with a subject scale or subscale
#' @param data     an \code{edsurvey.data.frame}, a \code{light.edsurvey.data.frame}, or an \code{edsurvey.data.frame.list}
#' @param groupA   an expression or character expression that defines a condition for the subset.
#'                 This subset will be compared to \code{groupB}. If not specified, it will define
#'                 a whole sample as in \code{data}.
#' @param groupB   an expression or character expression that defines a condition for the subset.
#'                 This subset will be compared to \code{groupA}. If not specified, it will define
#'                 a whole sample as in \code{data}. If set to \code{NULL}, estimates for the second group
#'                 will be dropped.
#' @param percentiles a numeric vector. The \code{gap} function calculates the
#'                    mean when this
#'                    argument is omitted or set to \code{NULL}. Otherwise,
#'                    the gap at the percentile given is calculated.
#' @param achievementLevel the achievement level(s) at which percentages
#'                         should be calculated
#' @param achievementDiscrete a logical indicating if the achievement level
#'                            specified in the \code{achievementLevel}
#'                            argument should be interpreted as discrete
#'                            so that
#'                            just the percentage in that particular achievement
#'                            level
#'                            will be included. Defaults to \code{FALSE}
#'                            so that
#'                            the percentage at or above that achievement level
#'                            will be
#'                            included in the percentage.
#' @param stDev a logical, set to \code{TRUE} to calculate the gap in standard deviations.
#' @param targetLevel a character string. When specified, calculates the gap in
#'                    the percentage of students at
#'                    \code{targetLevel} in the \code{variable} argument. This is useful for
#'                    comparing the gap in the percentage of students at a
#'                    survey response level.
#' @param weightVar a character indicating the weight variable to use.
#'                  See Details.
#' @param jrrIMax    a numeric value; when using the jackknife variance estimation method, the default estimation option, \code{jrrIMax=1}, uses the
#'                   sampling variance from the first plausible value as the component for sampling variance estimation. The \code{Vjrr}
#'                   term, or sampeling variance term, can be estimated with any number of plausible values, and values larger than the number of
#'                   plausible values on the survey (including \code{Inf}) will result in all plausible values being used.
#'                   Higher values of \code{jrrIMax} lead to longer computing times and more accurate variance estimates.
#' @param varMethod  deprecated parameter, \code{gap} always uses the jackknife variance estimation
#' @param dropOmittedLevels a logical value. When set to the default value of
#'                      \code{TRUE}, drops those levels of
#'                      all factor variables.
#'                      Use \code{print} on an \code{edsurvey.data.frame}
#'                      to see the omitted levels.
#' @param defaultConditions a logical value. When set to the default value
#'                          of \code{TRUE}, uses the default
#'                          conditions stored in \code{edsurvey.data.frame}
#'                          to subset the data.
#'                          Use \code{print} on an \code{edsurvey.data.frame}
#'                          to see the default conditions.
#' @param recode a list of lists to recode variables. Defaults to \code{NULL}.
#'               Can be set as
#'               \code{recode} \code{=} \code{list(var1} \code{=}
#'               \code{list(from} \code{=} \code{c("a",} \code{"b",}
#'               \code{"c"),} \code{to} \code{=} \code{"d"))}.
#' @param referenceDataIndex a numeric used only when the \code{data} argument is an
#'                           \code{edsurvey.data.frame.list},
#'                           indicating which dataset is the reference
#'                           dataset that other datasets are compared with.
#'                           Defaults to 1.
#' @param returnVarEstInputs a logical value; set to \code{TRUE} to return the
#'                           inputs to the jackknife and imputation variance
#'                           estimates which allows for the
#'                           computation
#'                           of covariances between estimates.
#' @param returnSimpleDoF    a logical value set to \code{TRUE} to return the degrees
#'                           of freedom for some statistics (see Value
#'                           section) that do not have a
#'                           \emph{t}-test; useful primarily for further computation
#' @param returnSimpleN      a logical value set to \code{TRUE} to add the count
#'                           (\emph{n}-size) of observations included in groups A and B
#'                           in the percentage object
#' @param returnNumberOfPSU a logical value set to \code{TRUE} to return the number of
#'                          PSUs used in the calculation
#' @param noCov set the covariances to zero in result
#' @param pctMethod a character that is one of \code{unbiased} or \code{simple}.
#'                  See the help for \code{\link{percentile}} for more information.
#' @param includeLinkingError a logical value set to \code{TRUE} to include the
#'                            linking error in variance estimation.
#'                            Standard errors (e.g., \code{diffAAse}, \code{diffBBse},
#'                            and \code{diffABABse}) and \emph{p}-values (e.g., \code{diffAApValue},
#'                            \code{diffBBpValue}, and \code{diffABABpValue}) would be adjusted for
#'                            comparisons between digitally based assessments (DBA) and
#'                            paper-based assessments (PBA) data.
#'                            This option is supported only for NAEP data.
#'
#' @param omittedLevels this argument is deprecated. Use \code{dropOmittedLevels}.
#'
#' @details This function calculates the gap between \code{groupA} and \code{groupB} (which
#' may be omitted to indicate the full sample). The gap is
#' calculated for one of four statistics:
#' \describe{
#'   \item{the gap in means}{The mean score gap (in the score
#'      variable) identified in the \code{variable} argument.
#'      This is the default. The means and their standard errors are
#'      calculated using the methods
#'      described in the \code{\link{lm.sdf}} function documentation.}
#'   \item{the gap in percentiles}{The gap between respondents at
#'      the percentiles specified in the \code{percentiles} argument.
#'      This is returned when the \code{percentiles} argument is
#'      defined. The mean and standard error are computed as described in the
#'      \code{\link{percentile}} function documentation.}
#'   \item{the gap in achievement levels}{The gap in the percentage of
#'      students at (when \code{achievementDiscrete} is \code{TRUE}) or at
#'      or above (when \code{achievementDiscrete} is \code{FALSE}) a
#'      particular achievement level. This is used when the
#'      \code{achievementLevel} argument is defined. The mean and standard error
#'      are calculated as described in the \code{\link{achievementLevels}}
#'      function documentation.}
#'   \item{the gap in a survey response}{The gap in the percentage of
#'      respondents responding at \code{targetLevel} to
#'      \code{variable}. This is used when \code{targetLevel} is
#'      defined. The mean and standard deviation are calculated as described in
#'      the \code{\link{edsurveyTable}} function documentation.}
#' }
#'
#' @return
#' The return type depends on if the class of the \code{data} argument is an
#' \code{edsurvey.data.frame} or an \code{edsurvey.data.frame.list}. Both
#' include the call (called \code{call}), a list called \code{labels},
#' an object named \code{percentage}
#' that shows the percentage in \code{groupA} and \code{groupB}, and an object
#' that shows the gap called \code{results}.
#'
#' The labels include the following elements:
#'   \item{definition}{the definitions of the groups}
#'   \item{nFullData}{the \emph{n}-size for the full dataset (before applying the definition)}
#'   \item{nUsed}{the \emph{n}-size for the data after the group is subsetted and other
#'                restrictions (such as omitted values) are applied}
#'   \item{nPSU}{the number of PSUs used in calculation--only returned when
#'               \code{returnNumberOfPSU} \code{=} \code{TRUE}}
#'
#' The percentages are computed according to the vignette titled
#' \href{https://www.air.org/sites/default/files/EdSurvey-Statistics.pdf}{\emph{Statistical Methods Used in EdSurvey}}
#'  in the section
#' \dQuote{Estimation of Weighted Percentages When Plausible Values Are Not Present.}
#' The standard errors are calculated according to
#' \dQuote{Estimation of the Standard Error of Weighted Percentages When Plausible Values Are Not Present, Using the Jackknife Method.}
#' Standard errors of differences are calculated as the square root of the typical
#' variance formula
#' \deqn{Var(A-B) = Var(A) + Var(B) - 2 Cov(A,B)}
#' where the covariance term is calculated as described in the vignette titled
#' \href{https://www.air.org/sites/default/files/EdSurvey-Statistics.pdf}{\emph{Statistical Methods Used in EdSurvey}}
#'  in the section
#' \dQuote{Estimation of Covariances.} These degrees of freedom are available only
#' with the jackknife variance estimation. The degrees of freedom used for hypothesis testing
#' are always set to the number of jackknife replicates in the data.
#'
#' \strong{the data argument is an edsurvey.data.frame}
#'   When the \code{data} argument is an \code{edsurvey.data.frame},
#'   \code{gap} returns an S3 object of class \code{gap}.
#'
#'   The \code{percentage} object is a numeric vector with the following elements:
#'     \item{pctA}{the percentage of respondents in \code{groupA} compared with the whole sample in \code{data}}
#'     \item{pctAse}{the standard error on the percentage of respondents in
#'                       \code{groupA}}
#'     \item{dofA}{degrees of freedom appropriate for a \emph{t}-test involving \code{pctA}.
#'                 This value is returned only if
#'                 \code{returnSimpleDoF}\code{=}\code{TRUE}.}
#'     \item{pctB}{the percentage of respondents in \code{groupB}.}
#'     \item{pctBse}{the standard error on the percentage of respondents in
#'                       \code{groupB}}
#'     \item{dofB}{degrees of freedom appropriate for a \emph{t}-test involving \code{pctA}.
#'                 This value is returned only if
#'                 \code{returnSimpleDoF}\code{=}\code{TRUE}.}
#'     \item{diffAB}{the value of \code{pctA} minus \code{pctB}}
#'     \item{covAB}{the covariance of \code{pctA} and \code{pctB}; used in
#'                  calculating \code{diffABse}.}
#'     \item{diffABse}{the standard error of \code{pctA}
#'                            minus \code{pctB}}
#'     \item{diffABpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffAB}
#'                         is zero}
#'     \item{dofAB}{degrees of freedom used in calculating
#'                       \code{diffABpValue}}
#'
#'   The \code{results} object is a numeric data frame with the following elements:
#'     \item{estimateA}{the mean estimate of \code{groupA} (or the percentage estimate
#'                      if \code{achievementLevel} or \code{targetLevel} is specified)}
#'     \item{estimateAse}{the standard error of \code{estimateA}}
#'     \item{dofA}{degrees of freedom appropriate for a \emph{t}-test involving \code{meanA}.
#'                 This value is returned only if
#'                 \code{returnSimpleDoF}\code{=}\code{TRUE}.}
#'     \item{estimateB}{the mean estimate of \code{groupB} (or the percentage estimate
#'                      if \code{achievementLevel} or \code{targetLevel} is specified)}
#'     \item{estimateBse}{the standard error of \code{estimateB}}
#'     \item{dofB}{degrees of freedom appropriate for a \emph{t}-test involving \code{meanB}.
#'                 This value is returned only if
#'                 \code{returnSimpleDoF}\code{=}\code{TRUE}.}
#'     \item{diffAB}{the value of \code{estimateA} minus \code{estimateB}}
#'     \item{covAB}{the covariance of \code{estimateA} and \code{estimateB}. Used in
#'                  calculating \code{diffABse}.}
#'     \item{diffABse}{the standard error of \code{diffAB}}
#'     \item{diffABpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffAB}
#'                         is zero.}
#'     \item{dofAB}{degrees of freedom used for the \emph{t}-test on \code{diffAB}}
#'
#'   If the gap was in  achievement levels or percentiles and more
#'   than one percentile or achievement level is requested,
#'   then an additional column
#'   labeled \code{percentiles} or \code{achievementLevel} is included
#'   in the \code{results} object.
#'
#'   When \code{results} has a single row and when \code{returnVarEstInputs}
#'   is \code{TRUE}, the additional elements \code{varEstInputs} and
#'   \code{pctVarEstInputs} also are returned. These can be used for calculating
#'   covariances with \code{\link{varEstToCov}}.
#'
#' \strong{the data argument is an edsurvey.data.frame.list}
#'   When the \code{data} argument is an \code{edsurvey.data.frame.list},
#'   \code{gap} returns an S3 object of class \code{gapList}.
#'
#'   The \code{results} object in the \code{edsurveyResultList} is
#'   a \code{data.frame}. Each row regards a particular dataset from the
#'   \code{edsurvey.data.frame}, and a reference dataset is dictated by
#'   the \code{referenceDataIndex} argument.
#'
#'   The \code{percentage} object is a \code{data.frame} with the following elements:
#'     \item{covs}{a data frame with a column for each column in the \code{covs}. See previous
#'                 section for more details.}
#'     \item{...}{all elements in the \code{percentage} object in the
#'                previous section}
#'     \item{diffAA}{the difference in \code{pctA} between the reference data
#'                   and this dataset. Set to \code{NA} for the
#'                   reference dataset.}
#'     \item{covAA}{the covariance of \code{pctA} in the reference data and
#'                  \code{pctA} on this row. Used in
#'                  calculating \code{diffAAse}.}
#'     \item{diffAAse}{the standard error for \code{diffAA}}
#'     \item{diffAApValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffAA}
#'                         is zero}
#'     \item{diffBB}{the difference in \code{pctB} between the reference data
#'                   and this dataset. Set to \code{NA} for the
#'                   reference dataset.}
#'     \item{covBB}{the covariance of \code{pctB} in the reference data and
#'                  \code{pctB} on this row. Used in
#'                  calculating \code{diffAAse}.}
#'     \item{diffBBse}{the standard error for \code{diffBB}}
#'     \item{diffBBpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffBB}
#'                         is zero}
#'     \item{diffABAB}{the value of \code{diffAB} in the reference dataset
#'                            minus the value of \code{diffAB} in this dataset. Set
#'                            to \code{NA} for the reference dataset.}
#'     \item{covABAB}{the covariance of \code{diffAB} in the reference data and
#'                    \code{diffAB} on this row. Used in
#'                    calculating \code{diffABABse}.}
#'     \item{diffABABse}{the standard error for \code{diffABAB}}
#'     \item{diffABABpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffABAB}
#'                         is zero}
#'
#'   The \code{results} object is a \code{data.frame} with the following elements:
#'     \item{...}{all elements in the \code{results} object in the
#'                previous section}
#'     \item{diffAA}{the value of \code{groupA} in the reference dataset minus
#'                          the value in this dataset. Set to \code{NA} for the
#'                          reference dataset.}
#'     \item{covAA}{the covariance of \code{meanA} in the reference data and
#'                  \code{meanA} on this row. Used in
#'                  calculating \code{diffAAse}.}
#'     \item{diffAAse}{the standard error for \code{diffAA}}
#'     \item{diffAApValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffAA}
#'                         is zero}
#'     \item{diffBB}{the value of \code{groupB} in the reference dataset minus
#'                          the value in this dataset. Set to \code{NA} for the
#'                          reference dataset.}
#'     \item{covBB}{the covariance of \code{meanB} in the reference data and
#'                  \code{meanB} on this row. Used in
#'                  calculating \code{diffBBse}.}
#'     \item{diffBBse}{the standard error for \code{diffBB}}
#'     \item{diffBBpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffBB}
#'                         is zero}
#'     \item{diffABAB}{the value of \code{diffAB} in the reference dataset
#'                            minus the value of \code{diffAB}
#'                            in this dataset. Set
#'                            to \code{NA} for the reference dataset.}
#'     \item{covABAB}{the covariance of \code{diffAB} in the reference data and
#'                    \code{diffAB} on this row. Used in
#'                    calculating \code{diffABABse}.}
#'     \item{diffABABse}{the standard error for \code{diffABAB}}
#'     \item{diffABABpValue}{the \emph{p}-value associated with the \emph{t}-test used
#'                         for the hypothesis test that \code{diffABAB}
#'                         is zero}
#'     \item{sameSurvey}{a logical value indicating if this line uses the same
#'                       survey as the reference line. Set to \code{NA} for the
#'                       reference line.}
#'
#' @author Paul Bailey, Trang Nguyen, and Huade Huo
#' @importFrom stats formula
#' @example man/examples/gap.R
#' @export
gap <- function(variable, data, groupA = "default", groupB = "default",
                percentiles = NULL,
                achievementLevel = NULL,
                achievementDiscrete = FALSE,
                stDev = FALSE,
                targetLevel = NULL,
                weightVar = NULL, jrrIMax = 1,
                varMethod = c("jackknife"),
                dropOmittedLevels = TRUE,
                defaultConditions = TRUE,
                recode = NULL,
                referenceDataIndex = 1,
                returnVarEstInputs = FALSE,
                returnSimpleDoF = FALSE,
                returnSimpleN = FALSE,
                returnNumberOfPSU = FALSE,
                noCov = FALSE,
                pctMethod = c("unbiased", "symmetric", "simple"),
                includeLinkingError = FALSE,
                omittedLevels = deprecated()) {
  if (is.character(substitute(groupA))) {
    groupA <- parse(text = groupA)[[1]]
  }
  if (all(as.character(substitute(groupA)) == "default")) {
    groupB <- NULL
  }
  if (!missing(groupB)) {
    if (is.character(substitute(groupB))) {
      groupB <- parse(text = groupB)[[1]]
    }
  }
  if (!missing(varMethod)) {
    warning(paste0("Argument ", dQuote("varMethod"), " deprecated. The ", dQuote("gap"), " function always uses jackknife variance estimation."))
  }

  if (lifecycle::is_present(omittedLevels)) {
    lifecycle::deprecate_soft("4.0.0", "gap(omittedLevels)", "gap(dropOmittedLevels)")
    dropOmittedLevels <- omittedLevels
  }
  if (!is.logical(dropOmittedLevels)) stop("The ", sQuote("dropOmittedLevels"), " argument must be logical.")

  # check incoming variables
  call <- match.call()
  # get valid pctMethod
  pctMethod <- match.arg(pctMethod)

  checkDataClass(data, c("edsurvey.data.frame", "light.edsurvey.data.frame", "edsurvey.data.frame.list"))

  # if recode list is provided, it will be added to userConditions at first (before subset in groupA and groupB)
  if (!is.null(recode)) {
    data <- recode.sdf(data, recode)
  }

  if (inherits(data, "edsurvey.data.frame") | inherits(data, "light.edsurvey.data.frame")) {
    survey <- getAttributes(data, "survey")
  } else {
    # edsurvey.data.frame.list
    robustSurvey <- function(dat) {
      res <- tryCatch(getAttributes(dat, "survey"),
        error = function(e) {
          NULL
        }
      )
      return(res)
    }
    survey <- unique(unlist(lapply(data$data, robustSurvey)))
  }

  # check if linking error is supported for this survey
  if (!missing(includeLinkingError) && includeLinkingError) {
    checkDataClass(data, c("edsurvey.data.frame.list"))
    if (!any(c("NAEP", "PISA") %in% survey)) {
      stop("the argument ", dQuote("includeLinkingError"), " can only be set to ", dQuote("TRUE"), " when EdSurvey suports linking error for the survey. Currently EdSurvey supports linking error for NAEP and PISA.")
    }
  }
  # Set dbaLabels as NULL for linking error with non ESDFL object or
  # includeLinkingError set as False
  dbaLabels <- NULL

  # gapStat required for PISA regardless of linking error
  if ("PISA" %in% survey) {
    if (!is.null(percentiles)) {
      gapStat <- "percentile"
    } else if (!is.null(achievementLevel)) {
      gapStat <- "AL"
    } else if (stDev) {
      gapStat <- "SD"
    } else {
      gapStat <- "Mean"
    }
  }
  # Get parameters for linking error
  if (!missing(includeLinkingError) && includeLinkingError) {
    # the following is only for NAEP which has linking error within survey
    if ("NAEP" %in% survey) {
      if (inherits(data, "edsurvey.data.frame.list")) {
        uniqueSubject <- unique(unlist(lapply(data$data, getAttributes, "subject")))
        if (length(uniqueSubject) != 1) {
          stop(paste0("Linking error only support ", dQuote("edsurvey.data.frame.list"), " with a single subject. Found subjects: ", pasteItems(uniqueSubject), "."))
        }

        uniqueGrade <- unique(unlist(lapply(data$data, getAttributes, "gradeLevel")))
        if (length(uniqueGrade) != 1) {
          stop(paste0("Linking error only support ", dQuote("edsurvey.data.frame.list"), " with a single grade level. Found grades: ", pasteItems(uniqueGrade), "."))
        }
      }

      # For achievementDiscrete=TRUE, return "At basic" (not "below basic") if "achievementLevel = "basic"
      alLinkingErrorhelper <- function(alstrings, discrete) {
        matchedAL <- NULL
        for (alstring in alstrings) {
          if (discrete) {
            if (grepl("^basic", alstring, ignore.case = TRUE)) {
              matchedAL <- c(matchedAL, "At basic")
            } else {
              matchedAL <- c(
                matchedAL,
                c(
                  "Below basic", "At basic", "At proficient",
                  "At advanced"
                )[grepl(alstring, c(
                  "Below basic", "At basic",
                  "At proficient", "At advanced"
                ),
                ignore.case = TRUE
                )]
              )
            }
          } else {
            matchedAL <- c(
              matchedAL,
              c(
                "At or above basic",
                "At or above proficient",
                "At advanced"
              )[
                grepl(alstring, c("At or above basic", "At or above proficient", "At advanced"),
                  ignore.case = TRUE
                )
              ]
            )
          }
        }
        if (length(matchedAL) != length(alstrings)) {
          stop("Linking error for provided achievement level ", pasteItems(alstrings[!alstrings %in% matchedAL]), "is not available.")
        }
        return(matchedAL)
      }

      # Get DBA labels from DBA years
      # 2017 is the first NAEP DBA year
      if (inherits(data, "edsurvey.data.frame.list")) {
        allYears <- unlist(itterateESDFL(call("getAttributes", attribute = "year"), data))
        dbaLabels <- data[[2]]$labels[allYears >= 2017]
        if (!all(dbaLabels %in% data[[2]]$labels)) {
          stop(paste0(
            pasteItems(dbaLabels[!dbaLabels %in% data[[2]]$labels]),
            " is not in the ", dQuote("edsurvey.data.frame.list"), " labels."
          ))
        }
      } else {
        allYears <- getAttributes(data, attribute = "year")
        dbaLabels <- NULL
      }
      if (!is.null(percentiles) && !all(percentiles %in% c(10, 25, 50, 75, 90))) {
        stop(
          "NAEP linking error only supported for 10th, 25th, 50th, 75th, 90th percentiles.",
          "Set ", dQuote("includeLinkingError"), " to ", dQuote("FALSE"),
          "to calculate any percentile without linking error."
        )
      }
      if (!is.null(percentiles)) {
        gapStat <- percentiles
      } else if (!is.null(achievementLevel)) {
        gapStat <- alLinkingErrorhelper(achievementLevel, achievementDiscrete)
      } else if (stDev) {
        gapStat <- "SD"
      } else {
        gapStat <- "Mean"
      }
    }
  }

  # deal with the possibility that there is more than one line of results for various possible reasons
  if (inherits(data, "edsurvey.data.frame.list") | length(percentiles) > 1 | length(achievementLevel) > 1) {
    # a flag that is TRUE when there is no information in data$covs. Default to false
    nocovs <- FALSE
    if (!inherits(data, "edsurvey.data.frame.list")) {
      # this is a situation where data is not an esdfl. The following code
      # required it to be one, so make it one. But note that covs
      # does not have real information in it.
      nocovs <- TRUE
      data <- edsurvey.data.frame.list(list(data), labels = "data")
    }
    ll <- length(data$datalist)
    lpct <- ifelse(is.null(percentiles), 1, length(percentiles))
    lal <- ifelse(is.null(achievementLevel), 1, length(achievementLevel))

    # check variable specific to edsurvey.data.frame.list
    suppressWarnings(refi <- as.integer(referenceDataIndex))
    if (!refi %in% 1:ll) {
      stop(paste0("The argument ", sQuote("referenceDataIndex"), " must be an integer between 1 and ", ll, "."))
    }
    res <- list(
      results = list(),
      labels = list(A = substitute(groupA), B = substitute(groupB))
    )

    pctdf0 <- data$covs
    if (!nocovs) {
      pctdf0$covAA <- 0
      pctdf0$dofAA <- Inf
      pctdf0$covBB <- 0
      pctdf0$dofBB <- Inf
      pctdf0$covABAB <- 0
      pctdf0$dofABAB <- Inf
    }

    # make sure refi gets calculated first
    llvec <- unique(c(refi, 1:ll))
    resilist <- list()
    # call gap for each edsurvey.data.frame
    for (i in llvec) {
      # for each element of the edsurvey.data.frame.list
      calli <- call
      # adjust the call per the needs of this call
      calli$data <- data$datalist[[i]]
      # we need the these inputs to calculate across essurvey.data.frame
      # covariances (but, only if they are the same dataset) an issue
      # that will be resolved later.
      calli$returnVarEstInputs <- TRUE
      # this returns the degrees of freedom for meanA and meanB
      # used for diffAA and diffBB, respectively
      calli$returnSimpleDoF <- TRUE
      # if these are NULL then setting them changes nothing
      calli$percentiles <- percentiles
      calli$achievementLevel <- achievementLevel
      calli <- as.list(calli)
      # use the same list of arguments but change the function to gapHelper
      calli[[1]] <- as.name("gapHelper")
      resilist[[i]] <- tryCatch(resi <- eval(as.call(calli)),
        error = function(cond) {
          message(paste0(
            "An error occurred while working on a dataset ",
            pasteItems(dQuote(data$covs[i, ])),
            ". The results from that dataset will be excluded. Error message: ",
            cond
          ))
          return(0)
        }
      )
      # if there was not an error
      if (!inherits(resilist[[i]], "gap")) {
        # there was an error
        if (i %in% referenceDataIndex) {
          stop("An error prevented processing of the reference dataset. This must be fixed before comparisons to the reference dataset can be made.")
        }
      } else {
        resi$call <- NULL

        # the following code is to clean up the results (adding some columns to make the code to extract results from gapHelper more consistent)
        # after this code, we should expect that each gap result will have:
        # 1. resi$results will have 'statisticsLevel' column: indicate percentile or achievementLevel, or 'data' if none is specified
        # 2. resi$varEstInputs: PV and JK data frame will have 'Level' column similar to 'statisticsLevel' column in resi$results
        if (is.null(percentiles) && is.null(achievementLevel)) {
          resi$results$statisticsLevel <- "data"
        }
        # if no or only 1 percentile/achievementLevel is specified
        if (length(percentiles) <= 1 && length(achievementLevel) <= 1) {
          resi$varEstInputs$JK$Level <- "data"
          if (!is.null(resi$varEstInputs$PV)) { # PV is null for type = "eq"
            resi$varEstInputs$PV$Level <- "data"
          }
        }
        # in gapHelper, type = 'pct' will return 'percentiles' column and type = 'al' will return 'achievementLevel' column
        # this replacement is to use 1 consistent name for the result extraction code below
        # it will be renamed back to original names after the extraction
        colnames(resi$results) <- gsub("achievementLevel|percentiles", "statisticsLevel", colnames(resi$results))
        resilist[[i]] <- resi
      }
    } # end for(i in llvec)

    # Extract results and compute AA or BB statistics
    lstats <- max(lpct, lal)
    for (j in 1:lstats) {
      resdf <- data$covs
      if (!nocovs) {
        # only generate these columns if there are multiple datasets
        resdf$covAA <- 0
        resdf$dofAA <- Inf
        resdf$covBB <- 0
        resdf$dofBB <- Inf
        resdf$covABAB <- 0
        resdf$dofABAB <- Inf
      }
      # to check whether the two samples are from the same survey or not
      # we need to use DOFCorrection and compute covariance if the two samples are from the same survey
      resdf$sameSurvey <- ifelse(1:nrow(resdf) == refi,
        rep(NA, nrow(resdf)),
        rep(FALSE, nrow(resdf))
      )
      if (lpct > 1) {
        resdf$percentiles <- percentiles[[j]]
      }
      for (i in llvec) {
        if (!inherits(resilist[[i]], "gap")) {
          # there was an error
          if (i %in% referenceDataIndex) {
            stop("An error prevented processing of the reference dataset. This must be fixed before comparisons to the reference dataset can be made.")
          }
          # this is to make sure that all achievementLevels or percentiles are created as to the reference edsurvey.data.frame
          # Missing values will be NA for statistics
          if (!is.null(achievementLevel)) {
            resdf$achievementLevel[i] <- resilist[[refi]]$results$statisticsLevel[j]
          } else if (!is.null(percentiles)) {
            resdf$percentiles[i] <- resilist[[refi]]$results$statisticsLevel[j]
          }
        } else { # there was not an error
          # dataframes do not like vectorwise assignment by column name,
          # so assign one at a time
          resi <- resilist[[i]]

          resi$results <- resi$results[j, ]
          resi$varEstInputs$JK <- resi$varEstInputs$JK[tolower(resi$varEstInputs$JK$Level) %in% tolower(c(resi$results$statisticsLevel, "data")), ]
          if (!is.null(resi$varEstInputs$PV)) {
            resi$varEstInputs$PV <- resi$varEstInputs$PV[tolower(resi$varEstInputs$PV$Level) %in% tolower(c(resi$results$statisticsLevel, "data")), ]
          }

          if (!is.null(achievementLevel)) {
            colnames(resi$results) <- gsub("statisticsLevel", "achievementLevel", colnames(resi$results))
          } else if (!is.null(percentiles)) {
            colnames(resi$results) <- gsub("statisticsLevel", "percentiles", colnames(resi$results))
          } else {
            resi$results$statisticsLevel <- NULL
          }
          for (ii in seq_along(resi$results)) {
            resdf[i, names(resi$results)[ii]] <- resi$results[[ii]]
          }
          skipB <- is.null(resi$labels$B) # a boolean to indicate whether groupB estimates are available in the results of gapHelper
          if (i == refi) {
            refVarEstInputs <- resi$varEstInputs
            refPctVarEstInputs <- resi$pctVarEstInputs
          }
          # refi is the index of reference dataset
          # AA and BB statistics of reference dataset will be set to NA
          if (i != refi) {
            resdf$dofAA[i] <- dofCompute((resi$results)$estimateAse, resdf[refi, "estimateAse"], (resi$results)[["dofA"]], resdf[refi, "dofA"])
            if (!skipB) {
              resdf$dofBB[i] <- dofCompute((resi$results)$estimateBse, resdf[refi, "estimateBse"], (resi$results)[["dofB"]], resdf[refi, "dofB"])
              resdf$dofABAB[i] <- dofCompute((resi$results)$diffABse, resdf[refi, "diffABse"], (resi$results)[["dofAB"]], resdf[refi, "dofAB"])
            }
            # if the datasets and from the same sample, add covariances and do DoFCorrection
            if (sameSurvey(data$datalist[[i]], data$datalist[[refi]])) {
              if (nrow(resi$varEstInputs$JK) != 0) { # some Level (percentile or achievementLevel) do not have varEstInput data so the subset
                # in line 476 will return an empty data.frame
                resdf$dofAA[i] <- DoFCorrection(resi$varEstInputs, refVarEstInputs, "A", method = "JR")
              }
              resdf$covAA[i] <- varEstToCov(resi$varEstInputs, refVarEstInputs, "A", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
              if (!skipB & nrow(resi$varEstInputs$JK) != 0) {
                resdf$dofBB[i] <- DoFCorrection(resi$varEstInputs, refVarEstInputs, "B", method = "JR")
                resdf$dofABAB[i] <- DoFCorrection(resi$varEstInputs, refVarEstInputs, "A-B", method = "JR")
                resdf$covBB[i] <- varEstToCov(resi$varEstInputs, refVarEstInputs, "B", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
                resdf$covABAB[i] <- varEstToCov(resi$varEstInputs, refVarEstInputs, "A-B", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
              }
              resdf$sameSurvey[i] <- TRUE
            } else {
              if ("PISA" %in% survey) {
                linke <- 0
                if (includeLinkingError & gapStat %in% c("Mean", "percentile", "SD")) {
                  if (variable %in% c("read", "math", "scie")) {
                    yearref <- getAttributes(data$datalist[[refi]], "year")
                    yeari <- getAttributes(data$datalist[[i]], "year")
                    if (yeari != yearref) {
                      linke <- calLinkingErrorPISA(subject = variable, years = c(yearref, yeari))
                    }
                  } else {
                    warning("No published linking error for variable ", dQuote(variable), "\n")
                  }
                }
                covvars <- c("covAA", "covBB", "covABAB")
                for (cvi in seq_along(covvars)) {
                  if (covvars[cvi] %in% colnames(resdf)) {
                    resdf[i, covvars[cvi]] <- -1 / 2 * linke^2
                  }
                }
              }
            }
          } else { # end if(i != refi)
            # if this is refi, set these to NULL
            resdf$covAA[i] <- resdf$covBB[i] <- resdf$covABAB[i] <- NA
            resdf$dofAA[i] <- resdf$dofBB[i] <- resdf$dofABAB[i] <- NA
          }
          # for the first j value only, grab the percent
          if (j == 1) {
            for (ii in 1:ncol(resi$percentage)) {
              pctdf0[i, colnames(resi$percentage)[ii]] <- resi$percentage[ , ii]
              if (i != refi) {
                pctdf0$dofAA[i] <- dofCompute((resi$percentage)$pctAse, pctdf0[refi, "pctAse"], (resi$percentage)[["dofA"]], pctdf0[refi, "dofA"])
                if (!skipB) {
                  pctdf0$dofBB[i] <- dofCompute((resi$percentage)$pctBse, pctdf0[refi, "pctBse"], (resi$percentage)[["dofB"]], pctdf0[refi, "dofB"])
                  pctdf0$dofABAB[i] <- dofCompute((resi$percentage)$diffABse, pctdf0[refi, "diffABse"], (resi$percentage)[["dofAB"]], pctdf0[refi, "dofAB"])
                }

                if (sameSurvey(data$datalist[[i]], data$datalist[[refi]])) {
                  pctdf0$dofAA[i] <- DoFCorrection(resi$pctVarEstInputs, refPctVarEstInputs, "A", method = "JR")
                  pctdf0$covAA[i] <- varEstToCov(resi$pctVarEstInputs, refPctVarEstInputs, "A", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
                  if (!skipB) {
                    pctdf0$dofBB[i] <- DoFCorrection(resi$pctVarEstInputs, refPctVarEstInputs, "B", method = "JR")
                    pctdf0$dofABAB[i] <- DoFCorrection(resi$pctVarEstInputs, refPctVarEstInputs, "A-B", method = "JR")
                    pctdf0$covBB[i] <- varEstToCov(resi$pctVarEstInputs, refPctVarEstInputs, "B", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
                    pctdf0$covABAB[i] <- varEstToCov(resi$pctVarEstInputs, refPctVarEstInputs, "A-B", jkSumMultiplier = getAttributes(data$datalist[[i]], "jkSumMultiplier"))
                  }
                }
              } else { # end if(i != refi)
                # if this is refi, set these to NULL
                pctdf0$covAA[i] <- pctdf0$covBB[i] <- pctdf0$covABAB[i] <- NA
                pctdf0$dofAA[i] <- pctdf0$dofBB[i] <- pctdf0$dofABAB[i] <- NA
              } # end else for (i != refi)
            } # end for(ii in 1:ncol(resi$percentage))
            if (returnSimpleN) {
              # return the n counts
              pctdf0$nA[i] <- resi$labels$nUsedA
              pctdf0$nB[i] <- resi$labels$nUsedB
            }
            if (returnNumberOfPSU) {
              pctdf0$nPSUA[i] <- ifelse(is.null(resi$labels$nPSUA), NA, resi$labels$nPSUA)
              pctdf0$nPSUB[i] <- ifelse(is.null(resi$labels$nPSUB), NA, resi$labels$nPSUB)
            }
          } # end if(j == 1)
        } # end if(inherits(temp, "gap"))
      } # End of for(i in 1:i) loop
      # Compute diff statistics (across different edsurvey.data.frame)
      resdf$diffAA <- ifelse(1:nrow(resdf) == refi, NA, resdf$estimateA[refi] - resdf$estimateA)
      if (!is.null(dbaLabels)) {
        # get linking error
        leAA <- calLinkingError(
          X = resdf$estimateA[refi],
          subjectI = uniqueSubject,
          gradeI = uniqueGrade,
          dependentVarI = variable,
          statElementI = gapStat,
          gapI = FALSE
        )[j]
        resdf$diffAAse <- ifelse(1:nrow(resdf) == refi, NA,
          ifelse(resdf$labels %in% dbaLabels,
            sqrt(resdf$estimateAse[refi]^2 + resdf$estimateAse^2 - 2 * resdf$covAA),
            sqrt(resdf$estimateAse[refi]^2 + resdf$estimateAse^2 - 2 * resdf$covAA + leAA)
          )
        )
      } else {
        resdf$diffAAse <- ifelse(1:nrow(resdf) == refi, NA, sqrt(resdf$estimateAse[refi]^2 + resdf$estimateAse^2 - 2 * resdf$covAA))
      }

      resdf$diffAApValue <- ifelse(1:nrow(resdf) == refi, NA, 2 * (1 - pt2(abs(resdf$diffAA / resdf$diffAAse), df = resdf$dofAA)))
      if (!skipB) {
        resdf$diffBB <- ifelse(1:nrow(resdf) == refi, NA, resdf$estimateB[refi] - resdf$estimateB)
        if (!is.null(dbaLabels)) {
          leBB <- calLinkingError(
            X = resdf$estimateB[refi],
            subjectI = uniqueSubject,
            gradeI = uniqueGrade,
            dependentVarI = variable,
            statElementI = gapStat,
            gapI = FALSE
          )[j]
          resdf$diffBBse <- ifelse(1:nrow(resdf) == refi, NA,
            ifelse(resdf$labels %in% dbaLabels,
              sqrt(resdf$estimateBse[refi]^2 + resdf$estimateBse^2 - 2 * resdf$covBB),
              sqrt(resdf$estimateBse[refi]^2 + resdf$estimateBse^2 - 2 * resdf$covBB + leBB)
            )
          )
        } else {
          resdf$diffBBse <- ifelse(1:nrow(resdf) == refi, NA, sqrt(resdf$estimateBse[refi]^2 + resdf$estimateBse^2 - 2 * resdf$covBB))
        }
        resdf$diffBBpValue <- ifelse(1:nrow(resdf) == refi, NA, 2 * (1 - pt2(abs(resdf$diffBB / resdf$diffBBse), df = resdf$dofBB)))
        resdf$diffABAB <- ifelse(1:nrow(resdf) == refi, NA, resdf$diffAB[refi] - resdf$diffAB)
        if (!is.null(dbaLabels)) {
          leABAB <- calLinkingError(
            X = resdf$diffAB[refi],
            subjectI = uniqueSubject,
            gradeI = uniqueGrade,
            dependentVarI = variable,
            statElementI = gapStat,
            gapI = TRUE
          )[j]
          resdf$diffABABse <- ifelse(1:nrow(resdf) == refi, NA,
            ifelse(resdf$labels %in% dbaLabels,
              sqrt(resdf$diffABse[refi]^2 + resdf$diffABse^2 - 2 * resdf$covABAB),
              sqrt(resdf$diffABse[refi]^2 + resdf$diffABse^2 - 2 * resdf$covABAB + leABAB)
            )
          )
        } else {
          resdf$diffABABse <- ifelse(1:nrow(resdf) == refi, NA, sqrt(resdf$diffABse[refi]^2 + resdf$diffABse^2 - 2 * resdf$covABAB))
        }
        resdf$diffABABpValue <- ifelse(1:nrow(resdf) == refi, NA, 2 * (1 - pt2(abs(resdf$diffABAB / resdf$diffABABse), df = resdf$dofABAB)))
      }

      # only generate these
      if (!nocovs) {
        pctdf0$diffAA <- ifelse(1:nrow(pctdf0) == refi, NA, pctdf0$pctA[refi] - pctdf0$pctA)
        pctdf0$diffAAse <- ifelse(1:nrow(pctdf0) == refi, NA, sqrt(pctdf0$pctAse[refi]^2 + pctdf0$pctAse^2 - 2 * pctdf0$covAA))
        pctdf0$diffAApValue <- ifelse(1:nrow(pctdf0) == refi, NA, 2 * (1 - pt2(abs(pctdf0$diffAA / pctdf0$diffAAse), df = pctdf0$dofAA)))
        if (!skipB) {
          pctdf0$diffBB <- ifelse(1:nrow(pctdf0) == refi, NA, pctdf0$pctB[refi] - pctdf0$pctB)
          pctdf0$diffBBse <- ifelse(1:nrow(pctdf0) == refi, NA, sqrt(pctdf0$pctBse[refi]^2 + pctdf0$pctBse^2 - 2 * pctdf0$covBB))
          pctdf0$diffBBpValue <- ifelse(1:nrow(pctdf0) == refi, NA, 2 * (1 - pt2(abs(pctdf0$diffBB / pctdf0$diffBBse), df = pctdf0$dofBB)))
          pctdf0$diffABAB <- ifelse(1:nrow(pctdf0) == refi, NA, pctdf0$diffAB[refi] - pctdf0$diffAB)
          pctdf0$diffABABse <- ifelse(1:nrow(pctdf0) == refi, NA, sqrt(pctdf0$diffABse[refi]^2 + pctdf0$diffABse^2 - 2 * pctdf0$covABAB))
          pctdf0$diffABABpValue <- ifelse(1:nrow(pctdf0) == refi, NA, 2 * (1 - pt2(abs(pctdf0$diffABAB / pctdf0$diffABABse), df = pctdf0$dofABAB)))
        }
      }
      if (j == 1) {
        resdf0 <- resdf
      } else {
        resdf0 <- rbind(resdf0, resdf)
      }
    } # end for(j in 1:lstats) loop
    if (nocovs) {
      # remove the labels columns which are just populated
      # with the unhelpful text "data"
      pctdf0$labels <- NULL
      resdf0$labels <- NULL
      # reorder columns when there are no covariates
      if (returnSimpleDoF) {
        if (!skipB) {
          resvar0 <- c(
            "estimateA", "estimateAse", "dofA",
            "estimateB", "estimateBse", "dofB"
          )
        } else {
          resvar0 <- c("estimateA", "estimateAse", "dofA")
        }
      } else {
        if (!skipB) {
          resvar0 <- c(
            "estimateA", "estimateAse",
            "estimateB", "estimateBse"
          )
        } else {
          resvar0 <- c("estimateA", "estimateAse")
        }
      }
      if (!skipB) {
        resvar0 <- c(
          resvar0,
          "diffAB", "covAB",
          "diffABse", "diffABpValue", "dofAB"
        )
      }
      addons <- c("achievementLevel", "percentiles")
      resvar <- c(addons[addons %in% names(resdf0)], resvar0)
      if ("sameSurvey" %in% names(resdf0) && sum(!is.na(resdf0$sameSurvey)) > 0) {
        resvar <- c(resvar, "sameSurvey")
      }
      resdf0 <- resdf0[ , resvar]
      # remove "estimate" with "pct" where "estimate" must start the variable name
      resvar <- gsub("^estimate", "pct", resvar0)
      if (returnSimpleN) {
        if (!skipB) {
          resvar <- c(resvar, "nA", "nB")
        } else {
          resvar <- c(resvar, "nA")
        }
      }
      if (returnNumberOfPSU) {
        if (!skipB) {
          resvar <- c(resvar, "nPSUA", "nPSUB")
        } else {
          resvar <- c(resvar, "nPSUA")
        }
      }
      pctdf0 <- pctdf0[ , resvar]
    } else { # end  if(nocovs)
      # reorder columns when there are covs
      if (returnSimpleDoF) {
        if (!skipB) {
          resvar0 <- c(
            names(data$covs), "estimateA", "estimateAse", "dofA",
            "estimateB", "estimateBse", "dofB"
          )
        } else {
          resvar0 <- c(names(data$covs), "estimateA", "estimateAse", "dofA")
        }
      } else {
        if (!skipB) {
          resvar0 <- c(
            names(data$covs), "estimateA", "estimateAse",
            "estimateB", "estimateBse"
          )
        } else {
          resvar0 <- c(names(data$covs), "estimateA", "estimateAse")
        }
      }
      if (!skipB) {
        resvar0 <- c(
          resvar0,
          "diffAB", "covAB", "diffABse",
          "diffABpValue", "dofAB",
          "diffAA", "covAA", "diffAAse",
          "diffAApValue", "dofAA",
          "diffBB", "covBB", "diffBBse",
          "diffBBpValue", "dofBB",
          "diffABAB", "covABAB", "diffABABse",
          "diffABABpValue", "dofABAB"
        )
      } else {
        resvar0 <- c(
          resvar0,
          "diffAA", "covAA", "diffAAse",
          "diffAApValue", "dofAA"
        )
      }
      addons <- c("achievementLevel", "percentiles")
      resvar <- c(addons[addons %in% names(resdf0)], resvar0)
      if ("sameSurvey" %in% names(resdf0) && sum(!is.na(resdf0$sameSurvey)) > 0) {
        resvar <- c(resvar, "sameSurvey")
      }
      resdf0 <- resdf0[ , resvar]
      resvar <- gsub("^estimate", "pct", resvar0)
      missing <- resvar[!resvar %in% names(pctdf0)]
      if (returnSimpleN) {
        if (!skipB) {
          resvar <- c(resvar, "nA", "nB")
        } else {
          resvar <- c(resvar, "nA")
        }
      }
      pctdf0 <- pctdf0[ , resvar]
    } # end else  if(nocovs)
    res[["results"]] <- resdf0
    res[["percentage"]] <- pctdf0
    res[["call"]] <- call
    class(res) <- "gapList"
    return(res)
  } else { # end if(inherits(data, "edsurvey.data.frame.list"))
    calli <- as.list(call)
    calli[[1]] <- as.name("gapHelper")
    calli$data <- data
    resi <- eval(as.call(calli))
    resi$call <- call
    return(resi)
  } # end else for if(inherits(data, "edsurvey.data.frame.list")) {
}


# this is called when data is edsurvey.data.frame or light.edsurvey.data.frame
gapHelper <- function(variable, data, groupA = "default", groupB = "default",
                      percentiles = NULL,
                      achievementLevel = NULL,
                      achievementDiscrete = FALSE,
                      stDev = FALSE,
                      targetLevel = NULL,
                      weightVar = NULL, jrrIMax = 1,
                      dropOmittedLevels = TRUE,
                      defaultConditions = TRUE,
                      recode = NULL,
                      referenceDataIndex = 1,
                      returnVarEstInputs = FALSE,
                      returnSimpleDoF = FALSE,
                      returnSimpleN = FALSE,
                      returnNumberOfPSU = FALSE,
                      noCov = FALSE,
                      pctMethod = c("unbiased", "symmetric", "simple"),
                      varMethod = NULL,
                      includeLinkingError = NULL) {
  if (is.character(substitute(groupA))) {
    groupA <- parse(text = groupA)[[1]]
  }
  if (all(as.character(substitute(groupA)) == "default")) {
    groupB <- NULL
  }
  if (!missing(groupB)) {
    if (is.character(substitute(groupB))) {
      groupB <- parse(text = groupB)[[1]]
    }
  }
  pctMethod <- match.arg(pctMethod)
  # get the weight var
  # if the weight var is not set, use the default
  wgt <- checkWeightVar(data, weightVar)

  varEstInputs <- NULL
  type <- "mu" # mean is the default
  if (!is.null(percentiles)) {
    type <- "pct"
  }
  if (!is.null(achievementLevel)) {
    if (type != "mu") {
      stop(paste0(
        "Only one of ", sQuote("percentiles"), ", ", sQuote("achievementLevel"),
        ", ", sQuote("SD"), ", and ", sQuote("targetLevel"), " can be defined."
      ))
    }
    type <- "AL"
    achievementLevel <- gmatchAttr(achievementLevel, getALNames(data, variable))
  }
  if (!is.null(targetLevel)) {
    if (type != "mu") {
      stop(paste0(
        "Only one of ", sQuote("percentiles"), ", ", sQuote("achievementLevel"),
        ", ", sQuote("SD"), ", and ", sQuote("targetLevel"), " can be defined."
      ))
    }
    type <- "eq"
  }
  if (stDev) {
    if (type != "mu") {
      stop(paste0(
        "Only one of ", sQuote("percentiles"), ", ", sQuote("achievementLevel"),
        ", ", sQuote("SD"), ", and ", sQuote("targetLevel"), " can be defined."
      ))
    }
    type <- "stdev"
  }
  # make the data with just groupA
  if (all(as.character(substitute(groupA)) == "default")) {
    dataA <- data
  } else {
    dataA <- subset(data, substitute(groupA), inside = TRUE)
  }
  # skipB is a boolean to indicate whether we need to compute statistics for group B
  # skipB = TRUE happens when (1) groupB is set to NULL or (2) groupA = 'default' (full sample)
  skipB <- FALSE
  # if groupB is null, just use data as the reference.
  # NOTE: the call to is.null requires "substitute" or groupB
  # will be evaluated right here
  if (is.null(substitute(groupB))) {
    skipB <- TRUE
  } else if (all(as.character(substitute(groupB)) == "default")) {
    dataB <- data
    dataAB <- dataA
  } else {
    # otherwise filter on groupB
    dataB <- subset(data, substitute(groupB), inside = TRUE)
    dataAB <- subset(dataA, substitute(groupB), inside = TRUE)
  }
  callv <- list(
    weightVar = wgt,
    jrrIMax = jrrIMax,
    dropOmittedLevels = dropOmittedLevels,
    recode = NULL
  ) # recode.sdf is already done in line 335
  if (!missing(defaultConditions)) {
    callv <- c(callv, list(defaultConditions = defaultConditions))
  }
  if (type == "mu") {
    # gap in mean
    calllm <- c(
      list(formula(paste0(variable, " ~ 1"))),
      callv,
      list(returnVarEstInputs = TRUE)
    ) # necessary for covariance est
    calllmA <- c(calllm, list(data = dataA))
    meanA <- do.call(lm.sdf, calllmA)
    coefA <- meanA$coefmat
    resA <- coefA[1, 1]
    if (!skipB) {
      calllmB <- c(calllm, list(data = dataB))
      meanB <- do.call(lm.sdf, calllmB)
      diff <- unname(meanA$coef[1] - meanB$coef[1]) # unname removes the name (Intercept) that otherwise ends up in the retruned value
      coefB <- meanB$coefmat
      resB <- coefB[1, 1]
      SEs <- list(estimateAse = coefA[1, 2], estimateBse = coefB[1, 2])
    } else {
      SEs <- list(estimateAse = coefA[1, 2])
    }

    # prepare varEstInputs
    # make varEstInputs$JK
    maJK <- subset(meanA$varEstInputs$JK, variable == "(Intercept)")
    maJK$variable <- "A"
    if (!skipB) {
      mbJK <- subset(meanB$varEstInputs$JK, variable == "(Intercept)")
      mbJK$variable <- "B"
      # calculate the difference
      mdJK <- merge(maJK, mbJK, by = c("PV", "JKreplicate"), suffixes = c(".A", ".B"))
      mdJK$variable <- "A-B"
      mdJK$value <- mdJK$value.A - mdJK$value.B
      mdJK$value.A <- mdJK$value.B <- NULL
      mdJK$variable.A <- mdJK$variable.B <- NULL
      varEstJK <- rbind(maJK, mbJK, mdJK)
    } else {
      varEstJK <- maJK
    }
    # make varEstInputs$PV
    if (!is.null(meanA$varEstInput$PV)) {
      maPV <- subset(meanA$varEstInputs$PV, variable == "(Intercept)")
      maPV$variable <- "A"
      if (!skipB) {
        mbPV <- subset(meanB$varEstInputs$PV, variable == "(Intercept)")
        mbPV$variable <- "B"
        # calculate the difference
        mdPV <- merge(maPV, mbPV, by = c("PV"), suffixes = c(".A", ".B"))
        mdPV$variable <- "A-B"
        mdPV$value <- mdPV$value.A - mdPV$value.B
        mdPV$value.A <- mdPV$value.B <- NULL
        mdPV$variable.A <- mdPV$variable.B <- NULL
        varEstPV <- rbind(maPV, mbPV, mdPV)
      } else {
        # only group A
        varEstPV <- maPV
      }
    } else { # end if(!is.null(meanA$varEstInput$PV))
      # no PVs, so no error from them
      varEstPV <- NULL
    } # end else for if(!is.null(meanA$varEstInput$PV))
    # make varEstInputs
    row.names(varEstJK) <- NULL
    row.names(varEstPV) <- NULL
    varEstInputs <- list(JK = varEstJK, PV = varEstPV)
  } # end if(type == "mu")
  if (type == "pct") {
    # gap at a percentile
    # here the length of percentiles is always one
    callpct <- c(
      list(
        variable = variable,
        percentiles = percentiles,
        confInt = FALSE, # do not share CI, so do not calculate it
        pctMethod = pctMethod
      ),
      callv,
      returnVarEstInputs = TRUE
    )
    callpctA <- c(callpct, list(data = dataA))
    meanA <- do.call(percentile, callpctA)
    resA <- meanA$estimate
    statisticsLevelOut <- paste0("P", percentiles)
    if (!skipB) {
      callpctB <- c(callpct, list(data = dataB))
      meanB <- do.call(percentile, callpctB)
      resB <- meanB$estimate
      SEs <- list(estimateAse = meanA$se, estimateBse = meanB$se)
    } else {
      SEs <- list(estimateAse = meanA$se)
    }
    maJK <- attributes(meanA)$varEstInputs$JK
    maJK$Level <- maJK$variable
    maJK$variable <- "A"
    if (!skipB) {
      mbJK <- attributes(meanB)$varEstInputs$JK
      mbJK$Level <- mbJK$variable
      mbJK$variable <- "B"
      # calculate the difference
      mdJK <- merge(maJK, mbJK, by = c("PV", "JKreplicate", "Level"), suffixes = c(".A", ".B"))
      mdJK$variable <- "A-B"
      mdJK$value <- mdJK$value.A - mdJK$value.B
      mdJK$value.A <- mdJK$value.B <- NULL
      mdJK$variable.A <- mdJK$variable.B <- NULL
      varEstJK <- rbind(maJK, mbJK, mdJK)
    } else {
      varEstJK <- maJK
    }

    # make varEstInputs$PV
    maPV <- attributes(meanA)$varEstInputs$PV
    maPV$Level <- maPV$variable
    maPV$variable <- "A"

    if (!skipB) {
      mbPV <- attributes(meanB)$varEstInputs$PV
      mbPV$Level <- mbPV$variable
      mbPV$variable <- "B"
      # calculate the difference
      mdPV <- merge(maPV, mbPV, by = c("PV", "Level"), suffixes = c(".A", ".B"))
      mdPV$variable <- "A-B"
      mdPV$value <- mdPV$value.A - mdPV$value.B
      mdPV$value.A <- mdPV$value.B <- NULL
      mdPV$variable.A <- mdPV$variable.B <- NULL
      varEstPV <- rbind(maPV, mbPV, mdPV)
    } else {
      varEstPV <- maPV
    }

    # make varEstInputs
    row.names(varEstJK) <- NULL
    row.names(varEstPV) <- NULL
    varEstInputs <- list(JK = varEstJK, PV = varEstPV)
  } # end if (type == 'pct')
  if (type == "AL") {
    # gap in percent at or above an achievement level
    callal <- c(
      list(
        achievementVars = variable,
        returnVarEstInputs = TRUE,
        returnCumulative = !achievementDiscrete
      ),
      callv
    )
    als <- getALNames(data, variable)
    callalA <- c(callal, list(data = dataA))
    meanA <- do.call(achievementLevels, callalA)
    if (!is.null(meanA$discrete)) {
      colnames(meanA$discrete)[grepl("Level$", colnames(meanA$discrete))] <- "Level"
    }
    if (!is.null(meanA$cumulative)) {
      colnames(meanA$cumulative)[grepl("Level$", colnames(meanA$cumulative))] <- "Level"
    }
    if (!skipB) {
      callalB <- c(callal, list(data = dataB))
      meanB <- do.call(achievementLevels, callalB)
      if (!is.null(meanB$discrete)) {
        colnames(meanB$discrete)[grepl("Level$", colnames(meanB$discrete))] <- "Level"
      }
      if (!is.null(meanB$cumulative)) {
        colnames(meanB$cumulative)[grepl("Level$", colnames(meanB$cumulative))] <- "Level"
      }
    }
    # achievementDiscrete indicates whether we should extract cumulative or discrete achievement level results
    if (achievementDiscrete) {
      # get index of desired achievement levels from the results
      if (typeof(als) == "character") {
        lA <- sapply(achievementLevel, function(al) {
          lal <- addALPrefix(al = al, als = als, discrete = TRUE)
          grep(lal, meanA$discrete$Level, ignore.case = TRUE)
        })
      } else {
        lA <- rep(TRUE, length(als[[achievementLevel]]))
      }

      if (length(lA) == 0) {
        stop("Achievement level cannot be found in the results of the call of the function ", sQuote("achievementLevel"), ".")
      }
      if (!skipB) {
        if (typeof(als) == "character") {
          lB <- sapply(achievementLevel, function(al) {
            lal <- addALPrefix(al = al, als = als, discrete = TRUE)
            grep(lal, meanB$discrete$Level, ignore.case = TRUE)
          })
        } else {
          lB <- rep(TRUE, length(als[[achievementLevel]]))
        }
        if (length(lB) == 0) {
          stop("Achievement level cannot be found in the results of the call of the function ", sQuote("achievementLevel"), ".")
        }
      }
      resA <- meanA$discrete$Percent[unlist(lA)]
      if (!skipB) {
        resB <- meanB$discrete$Percent[unlist(lB)]
        SEs <- list(
          estimateAse = meanA$discrete$StandardError[unlist(lA)],
          estimateBse = meanB$discrete$StandardError[unlist(lB)]
        )
      } else {
        SEs <- list(estimateAse = meanA$discrete$StandardError[unlist(lA)])
      }
    } else {
      # get index of desired achievement levels from the results
      if (typeof(als) == "character") {
        lA <- sapply(achievementLevel, function(al) {
          lal <- addALPrefix(al = al, als = als, discrete = FALSE)
          grep(lal, meanA$cumulative$Level, ignore.case = TRUE)
        })
      } else {
        lA <- rep(TRUE, length(als[[achievementLevel]]))
      }
      if (length(lA) == 0) {
        stop("Achievement level cannot be found in the results of the call of the function ", sQuote("achievementLevel"), ".")
      }
      if (!skipB) {
        if (typeof(als) == "character") {
          lB <- sapply(achievementLevel, function(al) {
            lal <- addALPrefix(al = al, als = als, discrete = FALSE)
            grep(lal, meanB$cumulative$Level, ignore.case = TRUE)
          })
        } else {
          lB <- rep(TRUE, length(als[[achievementLevel]]))
        }
      }
      resA <- meanA$cumulative$Percent[unlist(lA)]
      if (!skipB) {
        resB <- meanB$cumulative$Percent[unlist(lB)]
        SEs <- list(
          estimateAse = meanA$cumulative$StandardError[unlist(lA)],
          estimateBse = meanB$cumulative$StandardError[unlist(lB)]
        )
      } else {
        SEs <- list(estimateAse = meanA$cumulative$StandardError[unlist(lA)])
      }
    }
    # used inside next if, but also outside
    statisticsLevelOut <- meanA$discrete$Level[unlist(lA)]
    # get Cov
    if (achievementDiscrete) {
      maJK <- subset(meanA$discVarEstInputs$JK, variable %in% meanA$discrete$Level[unlist(lA)])
      maJK$Level <- maJK$variable
      maJK$variable <- "A"
      if (!skipB) {
        mbJK <- subset(meanB$discVarEstInputs$JK, variable %in% meanB$discrete$Level[unlist(lB)])
        mbJK$Level <- mbJK$variable
        mbJK$variable <- "B"
      }
    } # end if(achievementDiscrete)
    else {
      maJK <- subset(meanA$cumVarEstInputs$JK, variable %in% meanA$cumulative$Level[unlist(lA)])
      maJK$Level <- maJK$variable
      maJK$variable <- "A"
      if (!skipB) {
        mbJK <- subset(meanB$cumVarEstInputs$JK, variable %in% meanB$cumulative$Level[unlist(lB)])
        mbJK$Level <- mbJK$variable
        mbJK$variable <- "B"
      }
    } # end else if(achievementDiscrete)
    if (!skipB) {
      mdJK <- merge(maJK, mbJK, by = c("PV", "JKreplicate", "Level"), suffixes = c(".A", ".B"))
      mdJK$variable <- "A-B"
      mdJK$value <- mdJK$value.A - mdJK$value.B
      mdJK$value.A <- mdJK$value.B <- NULL
      mdJK$variable.A <- mdJK$variable.B <- NULL
      varEstJK <- rbind(maJK, mbJK, mdJK)
    } else {
      varEstJK <- maJK
    }
    # make varEstInputs$PV
    if (achievementDiscrete) {
      maPV <- subset(meanA$discVarEstInputs$PV, variable %in% meanA$discrete$Level[unlist(lA)])
      maPV$Level <- maPV$variable
      maPV$variable <- "A"
      if (!skipB) {
        mbPV <- subset(meanB$discVarEstInputs$PV, variable %in% meanB$discrete$Level[unlist(lB)])
        mbPV$Level <- mbPV$variable
        mbPV$variable <- "B"
      }
    } # end if(achievementDiscrete)
    else {
      statisticsLevelOut <- meanA$cumulative$Level[unlist(lA)]
      maPV <- subset(meanA$cumVarEstInputs$PV, variable %in% meanA$cumulative$Level[unlist(lA)])
      maPV$Level <- maPV$variable
      maPV$variable <- "A"
      if (!skipB) {
        mbPV <- subset(meanB$cumVarEstInputs$PV, variable %in% meanB$cumulative$Level[unlist(lB)])
        mbPV$Level <- mbPV$variable
        mbPV$variable <- "B"
      }
    } # end else if(achievementDiscrete)
    # calculate the difference
    if (!skipB) {
      mdPV <- merge(maPV, mbPV, by = c("PV", "Level"), suffixes = c(".A", ".B"))
      mdPV$variable <- "A-B"
      mdPV$value <- mdPV$value.A - mdPV$value.B
      mdPV$value.A <- mdPV$value.B <- NULL
      mdPV$variable.A <- mdPV$variable.B <- NULL
      varEstPV <- rbind(maPV, mbPV, mdPV)
    } else {
      varEstPV <- maPV
    }
    # make varEstInputs
    row.names(varEstJK) <- NULL
    row.names(varEstPV) <- NULL
    varEstInputs <- list(JK = varEstJK, PV = varEstPV)
  } # end if (type == "AL")
  if (type == "stdev") {
    # gap at a standard deviation
    # here the length of standard deviation is always one
    callStDev <- c(list(variable = variable),
      callv,
      returnVarEstInputs = TRUE
    )
    if ("weightVar" %in% names(callStDev) && is.null(callStDev[["weightVar"]])) {
      # default is used only if weightVar is missing.
      callStDev[["weightVar"]] <- NULL
    }
    callStDevA <- c(callStDev, list(data = dataA))
    meanA <- do.call(SD, callStDevA)
    resA <- meanA$std
    statisticsLevelOut <- "SD"
    if (!skipB) {
      callStDevB <- c(callStDev, list(data = dataB))
      meanB <- do.call(SD, callStDevB)
      resB <- meanB$std
      SEs <- list(estimateAse = meanA$stdSE, estimateBse = meanB$stdSE)
    } else {
      SEs <- list(estimateAse = meanA$stdSE)
    }
    maJK <- meanA$varEstInputs$JK
    maJK$Level <- maJK$variable
    maJK$variable <- "A"
    if (!skipB) {
      mbJK <- meanB$varEstInputs$JK
      mbJK$Level <- mbJK$variable
      mbJK$variable <- "B"
      # calculate the difference
      mdJK <- merge(maJK, mbJK, by = c("PV", "JKreplicate", "Level"), suffixes = c(".A", ".B"))
      mdJK$variable <- "A-B"
      mdJK$value <- mdJK$value.A - mdJK$value.B
      mdJK$value.A <- mdJK$value.B <- NULL
      mdJK$variable.A <- mdJK$variable.B <- NULL
      varEstJK <- rbind(maJK, mbJK, mdJK)
    } else {
      varEstJK <- maJK
    }

    # make varEstInputs$PV
    if (!is.null(meanA$varEstInput$PV)) {
      maPV <- meanA$varEstInputs$PV
      maPV$Level <- maPV$variable
      maPV$variable <- "A"

      if (!skipB) {
        mbPV <- meanB$varEstInputs$PV
        mbPV$Level <- mbPV$variable
        mbPV$variable <- "B"
        # calculate the difference
        mdPV <- merge(maPV, mbPV, by = c("PV", "Level"), suffixes = c(".A", ".B"))
        mdPV$variable <- "A-B"
        mdPV$value <- mdPV$value.A - mdPV$value.B
        mdPV$value.A <- mdPV$value.B <- NULL
        mdPV$variable.A <- mdPV$variable.B <- NULL
        varEstPV <- rbind(maPV, mbPV, mdPV)
      } else {
        varEstPV <- maPV
      }
    } else { # end if(!is.null(meanA$varEstInput$PV))
      # no PVs, so no error from them
      varEstPV <- NULL
    } # end else for if(!is.null(meanA$varEstInput$PV))

    # make varEstInputs
    row.names(varEstJK) <- NULL
    row.names(varEstPV) <- NULL
    varEstInputs <- list(JK = varEstJK, PV = varEstPV)
  } # end if (type == 'stDev')
  groupA_0 <- FALSE
  groupB_0 <- FALSE
  if (type == "eq") {
    # gap in percent in some covariate
    calleq <- c(callv, list(
      formula = formula(paste0(" ~ ", variable)),
      returnMeans = FALSE,
      pctAggregationLevel = NULL,
      returnSepct = TRUE,
      drop = FALSE,
      returnVarEstInputs = TRUE
    ))
    calleqA <- c(calleq, list(data = dataA))
    meanA <- do.call(edsurveyTable, calleqA)

    if (!targetLevel %in% meanA$data[[variable]]) {
      groupA_0 <- TRUE
    }
    resA <- ifelse(groupA_0, 0, meanA$data[meanA$data[variable] == targetLevel, 4])
    if (!skipB) {
      calleqB <- c(calleq, list(data = dataB))
      meanB <- do.call(edsurveyTable, calleqB)
      if (!targetLevel %in% meanB$data[[variable]]) {
        groupB_0 <- TRUE
      }
      resB <- ifelse(groupB_0, 0, meanB$data[meanB$data[variable] == targetLevel, 4])
      SEs <- list(
        estimateAse = ifelse(groupB_0, NA, meanA$data[meanA$data[variable] == targetLevel, 5]),
        estimateBse = ifelse(groupB_0, NA, meanB$data[meanB$data[variable] == targetLevel, 5])
      )
    } else {
      SEs <- list(estimateAse = ifelse(groupB_0, NA, meanA$data[meanA$data[variable] == targetLevel, 5]))
    }

    # make varEstInputs$JK
    if (!groupA_0) {
      maJK <- meanA$pctVarEstInputs$JK[meanA$pctVarEstInputs$JK$variable == paste0(variable, "=", targetLevel), ]
      maJK$variable <- "A"
    } else {
      maJK <- NULL
    }
    if (!skipB) {
      if (!groupB_0) {
        mbJK <- meanB$pctVarEstInputs$JK[meanB$pctVarEstInputs$JK$variable == paste0(variable, "=", targetLevel), ]
        mbJK$variable <- "B"
      } else {
        mbJK <- NULL
      }
      # calculate the difference
      if (!groupA_0 & !groupB_0) {
        mdJK <- merge(maJK, mbJK, by = c("PV", "JKreplicate"), suffixes = c(".A", ".B"))
        mdJK$variable <- "A-B"
        mdJK$value <- mdJK$value.A - mdJK$value.B
        mdJK$value.A <- mdJK$value.B <- NULL
        mdJK$variable.A <- mdJK$variable.B <- NULL
      } else {
        mdJK <- NULL
      }
      varEstJK <- rbind(maJK, mbJK, mdJK)
    } else {
      varEstJK <- maJK
    }

    # make varEstInputs
    if (!is.null(varEstJK)) {
      row.names(varEstJK) <- NULL
      varEstInputs <- list(JK = varEstJK, PV = NULL)
    } else {
      varEstInputs <- NULL
    }
  } # end if (type == "eq")
  # Compute difference statistics (diffAB, covAB, pooledse, pooleddf)
  if (!skipB) {
    diff <- unname(resA - resB) # unname removes the name (Intercept) that otherwise ends up in the retruned value
    nSize <- c(nrow(dataA), nrow(dataB))
    total <- which.max(nSize) # total is 1 if groupA is the total and 2 if groupB is the total
    p <- nrow(dataAB) / nSize[total]
    cov <- p * (SEs[[total]])^2
    wgtl <- getAttributes(data, "weights")[[wgt]]
    JRdfA <- JRdfB <- JRdf <- length(wgtl$jksuffixes)
    if (!is.null(varEstInputs)) {
      if (type %in% c("pct", "AL")) {
        cov <- c()
        JRdfA <- c()
        JRdfB <- c()
        JRdf <- c()
        # we need to calculate cov by each Level
        for (si in statisticsLevelOut) {
          varEstInputsTemp <- list(
            JK = subset(varEstInputs$JK, Level %in% si),
            PV = subset(varEstInputs$PV, Level %in% si)
          )
          if (nrow(varEstInputsTemp$JK) == 0) {
            varEstInputsTemp$JK <- NULL
          }
          if (nrow(varEstInputsTemp$PV) == 0) {
            varEstInputsTemp$PV <- NULL
          }
          if (!is.null(varEstInputsTemp$JK)) {
            cov <- c(cov, varEstToCov(varEstInputsTemp,
              varA = "A", varB = "B", jkSumMultiplier = getAttributes(data, "jkSumMultiplier")
            ))

            JRdfA <- c(JRdfA, DoFCorrection(varEstInputsTemp, varA = "A", method = c("JR")))
            JRdfB <- c(JRdfB, DoFCorrection(varEstInputsTemp, varA = "B", method = c("JR")))
            JRdf <- c(JRdf, DoFCorrection(varEstInputsTemp, varA = "A", varB = "B", method = c("JR")))
          } # end if (!is.null(varEstInputsTemp$JK))
          else {
            cov <- c(cov, NA)
            JRdfA <- c(JRdfA, NA)
            JRdfB <- c(JRdfB, NA)
            JRdf <- c(JRdf, NA)
          }
        }
      } # end if (type %in% c("pct","AL"))
      else {
        cov <- ifelse(groupA_0 | groupB_0, cov, varEstToCov(varEstInputs, varA = "A", varB = "B", jkSumMultiplier = getAttributes(data, "jkSumMultiplier")))
        JRdfA <- ifelse(groupA_0, NA, DoFCorrection(varEstInputs, varA = "A", method = c("JR")))
        JRdfB <- ifelse(groupB_0, NA, DoFCorrection(varEstInputs, varA = "B", method = c("JR")))
        JRdf <- ifelse(groupA_0 | groupB_0, NA, DoFCorrection(varEstInputs, varA = "A", varB = "B", method = c("JR")))
      }
    } # end if (!is.null(varEstInputs))
    if (noCov) {
      cov <- 0 * cov
      JRdf <- ifelse(groupA_0 | groupB_0, NA, (SEs[[1]]^2 + SEs[[2]]^2)^2 / (SEs[[1]]^4 / JRdfA + SEs[[2]]^4 / JRdfB))
    }
    diffSE <- mapply(function(a, b, c) {
      unname(sqrt(a^2 + b^2 - 2 * c))
    }, SEs[[1]], SEs[[2]], cov)
    diffSE <- unlist(diffSE)
    # get percent in each group

    # get variables used in conditions for groupA and groupB
    vn <- c(parseVars(substitute(groupA), data), parseVars(substitute(groupB), data))
    # add weights and variable itself
    vn <- c(vn, wgt, variable)

    tryCatch(vn <- unique(c(vn, PSUStratumNeeded(returnNumberOfPSU, data))),
      error = function(e) {
        warning(paste0("Stratum and PSU variables are required for this call and are not on the incoming data. Ignoring ", dQuote("returnNumberOfPSU=TRUE"), "."))
        returnNumberOfPSU <<- FALSE
      }
    )

    # setup non-data call variables
    callv <- list(
      varnames = vn,
      drop = FALSE,
      dropOmittedLevels = dropOmittedLevels,
      recode = NULL
    ) # recode is already done in line 335
    if (!missing(defaultConditions)) {
      callv <- c(callv, list(defaultConditions = defaultConditions))
    }

    # add weight call variables and do calls
    wgtcall0 <- c(callv, list(data = data))
    d0 <- do.call(getData, wgtcall0)
    wgtcallA <- c(callv, list(data = dataA))
    dA <- do.call(getData, wgtcallA)
    wgtcallB <- c(callv, list(data = dataB))
    dB <- do.call(getData, wgtcallB)
    pctA <- sum(dA[ , wgt]) / sum(d0[ , wgt])
    pctB <- sum(dB[ , wgt]) / sum(d0[ , wgt])
    # calculate covAB using JK method
    wgtl <- getAttributes(data, "weights")[[wgt]]
    pctJK <- t(sapply(wgtl$jksuffixes, function(suffix) {
      w0 <- sum(d0[ , paste0(wgtl$jkbase, suffix)])
      wA <- sum(dA[ , paste0(wgtl$jkbase, suffix)])
      wB <- sum(dB[ , paste0(wgtl$jkbase, suffix)])
      c(pctA = wA / w0, pctB = wB / w0)
    }))
    seA <- 100 * sqrt(getAttributes(data, "jkSumMultiplier") * sum((pctJK[ , 1] - pctA)^2))
    seB <- 100 * sqrt(getAttributes(data, "jkSumMultiplier") * sum((pctJK[ , 2] - pctB)^2))
    varJK <- data.frame(
      stringsAsFactors = FALSE,
      PV = rep(0, 3 * nrow(pctJK)),
      JKreplicate = rep(1:nrow(pctJK), 3),
      variable = rep(c("A", "B", "A-B"), each = nrow(pctJK)),
      value = c(
        pctJK[ , 1] - pctA,
        pctJK[ , 2] - pctB,
        (pctJK[ , 1] - pctA) - (pctJK[ , 2] - pctB)
      )
    )
    # there is no PV here because weights are not imputed
    pctVarEstInputs <- list(JK = varJK, PV = NULL)
    if (noCov) {
      covAB <- rep(0, length(seA))
      pctDoFA <- DoFCorrection(pctVarEstInputs, varA = "A", method = "JR")
      pctDoFB <- DoFCorrection(pctVarEstInputs, varA = "B", method = "JR")
      # Satterthwaite approximation
      pctJRdf <- (seA^2 + seB^2)^2 / (seA^4 / pctDoFA + seB^4 / pctDoFB)
    } else {
      covAB <- varEstToCov(pctVarEstInputs, varA = "A", varB = "B", jkSumMultiplier = getAttributes(data, "jkSumMultiplier"))
      pctDoFA <- DoFCorrection(pctVarEstInputs, varA = "A", method = "JR")
      pctDoFB <- DoFCorrection(pctVarEstInputs, varA = "B", method = "JR")
      pctJRdf <- DoFCorrection(pctVarEstInputs, varA = "A-B", method = "JR")
    }
    # end covAB calculation
    pctA <- 100 * pctA
    pctB <- 100 * pctB
    seAB <- unname(sqrt(seA^2 + seB^2 - 2 * 100 * 100 * covAB))

    # start return vector, first name achievement level
    vec <- switch(type,
      pct = data.frame(stringsAsFactors = FALSE, percentiles = as.numeric(percentiles)),
      AL = data.frame(stringsAsFactors = FALSE, achievementLevel = gsub("Below", "below", as.character(statisticsLevelOut))),
      data.frame()
    )
    if (returnSimpleDoF) {
      # add dofA and dofB if requested
      if (nrow(vec) > 0) {
        vec <- cbind(
          vec,
          data.frame(
            stringsAsFactors = FALSE,
            estimateA = resA, estimateAse = SEs[[1]], dofA = JRdfA,
            estimateB = resB, estimateBse = SEs[[2]], dofB = JRdfB
          )
        )
      } else { # this is equivalent to  if (type %in%  c("mu","eq"))
        vec <- data.frame(
          stringsAsFactors = FALSE,
          estimateA = resA, estimateAse = SEs[[1]], dofA = JRdfA,
          estimateB = resB, estimateBse = SEs[[2]], dofB = JRdfB
        )
      }
    } # end if(returnSimpleDoF)
    else {
      # otherwise, no dofA or dofB
      if (nrow(vec) > 0) {
        vec <- cbind(
          vec,
          data.frame(
            stringsAsFactors = FALSE,
            estimateA = resA, estimateAse = SEs[[1]],
            estimateB = resB, estimateBse = SEs[[2]]
          )
        )
      } else { # this is equivalent to  if (type %in%  c("mu","eq"))
        vec <- data.frame(
          stringsAsFactors = FALSE,
          estimateA = resA, estimateAse = SEs[[1]],
          estimateB = resB, estimateBse = SEs[[2]]
        )
      }
    }
    diffp <- 2 * (1 - pt2(abs(diff / diffSE), df = JRdf))
    vec <- cbind(
      vec,
      data.frame(
        stringsAsFactors = FALSE,
        diffAB = diff, covAB = unname(cov),
        diffABse = diffSE,
        diffABpValue = diffp,
        dofAB = JRdf
      )
    )
    row.names(vec) <- NULL
    if (returnSimpleDoF) {
      percentage <- data.frame(
        stringsAsFactors = FALSE,
        pctA = pctA,
        pctAse = seA,
        dofA = pctDoFA,
        pctB = pctB,
        pctBse = seB,
        dofB = pctDoFB
      )
    } else {
      percentage <- data.frame(
        stringsAsFactors = FALSE,
        pctA = pctA,
        pctAse = seA,
        pctB = pctB,
        pctBse = seB
      )
    }
    pdiffp <- 2 * (1 - pt2(abs(pctA - pctB) / seAB, df = pctJRdf))
    percentage <- cbind(
      percentage,
      data.frame(
        stringsAsFactors = FALSE,
        diffAB = pctA - pctB,
        covAB = 100 * 100 * covAB,
        diffABse = seAB,
        diffABpValue = pdiffp,
        dofAB = pctJRdf
      )
    )

    lst <- list(
      results = vec,
      labels = list(
        A = substitute(groupA), B = substitute(groupB),
        n0A = nrow2.edsurvey.data.frame(dataA),
        n0B = nrow2.edsurvey.data.frame(dataB),
        nUsedA = nrow(dA),
        nUsedB = nrow(dB)
      ),
      percentage = percentage
    )
    if (returnNumberOfPSU) {
      lst$labels$nPSUA <- nrow(unique(dA[ , c(getAttributes(dataA, "stratumVar"), getAttributes(dataA, "psuVar"))]))
      lst$labels$nPSUB <- nrow(unique(dB[ , c(getAttributes(dataB, "stratumVar"), getAttributes(dataB, "psuVar"))]))
    }
  } else { # end if(!skipB)
    wgtl <- getAttributes(data, "weights")[[wgt]]
    JRdfA <- ifelse(groupA_0, NA, length(wgtl$jksuffixes))
    if (!is.null(varEstInputs)) {
      if (type %in% c("pct", "AL")) {
        JRdfA <- c()
        for (si in statisticsLevelOut) {
          varEstInputsTemp <- list(
            JK = subset(varEstInputs$JK, Level %in% si),
            PV = subset(varEstInputs$PV, Level %in% si)
          )
          if (nrow(varEstInputsTemp$JK) == 0) {
            varEstInputsTemp$JK <- NULL
          }
          if (nrow(varEstInputsTemp$PV) == 0) {
            varEstInputsTemp$PV <- NULL
          }
          if (!is.null(varEstInputsTemp$JK)) {
            JRdfA <- c(JRdfA, DoFCorrection(varEstInputsTemp, varA = "A", method = c("JR")))
          } else {
            JRdfA <- c(JRdfA, NA)
          }
        }
      } # end if (type %in% c("pct","AL"))
      else {
        JRdfA <- ifelse(groupA_0, NA, DoFCorrection(varEstInputs, varA = "A", method = c("JR")))
      }
    }
    # get variables used in conditions for groupA and groupB
    vn <- c(parseVars(substitute(groupA), data), parseVars(substitute(groupB), data))
    # add weights and variable itself
    vn <- c(vn, wgt, variable)
    tryCatch(vn <- unique(c(vn, PSUStratumNeeded(returnNumberOfPSU, data))),
      error = function(e) {
        warning(paste0("Stratum and PSU variables are required for this call and are not on the incoming data. Ignoring ", dQuote("returnNumberOfPSU=TRUE"), "."))
        returnNumberOfPSU <<- FALSE
      }
    )
    # setup non-data call variables
    callv <- list(
      varnames = vn,
      drop = FALSE,
      dropOmittedLevels = dropOmittedLevels,
      recode = NULL,
      addAttributes = TRUE
    ) # recode is already done in line 335
    if (!missing(defaultConditions)) {
      callv <- c(callv, list(defaultConditions = defaultConditions))
    }

    # add weight call variables and do calls
    wgtcall0 <- c(callv, list(data = data))
    d0 <- do.call(getData, wgtcall0)
    wgtcallA <- c(callv, list(data = dataA))
    dA <- do.call(getData, wgtcallA)
    pctA <- sum(dA[ , wgt]) / sum(d0[ , wgt])
    pctJK <- t(sapply(wgtl$jksuffixes, function(suffix) {
      w0 <- sum(d0[ , paste0(wgtl$jkbase, suffix)])
      wA <- sum(dA[ , paste0(wgtl$jkbase, suffix)])
      c(pctA = wA / w0, pctB = 0)
    }))
    seA <- 100 * sqrt(getAttributes(data, "jkSumMultiplier") * sum((pctJK[ , 1] - pctA)^2))
    varJK <- data.frame(
      stringsAsFactors = FALSE,
      PV = rep(0, nrow(pctJK)),
      JKreplicate = 1:nrow(pctJK),
      variable = "A",
      value = pctJK[ , 1] - pctA
    )
    # there is no PV here because weights are not imputed
    pctVarEstInputs <- list(JK = varJK, PV = NULL)
    pctDoFA <- DoFCorrection(pctVarEstInputs, varA = "A", method = "JR")
    pctA <- 100 * pctA
    # start return vector, first name achievement level
    vec <- switch(type,
      pct = data.frame(stringsAsFactors = FALSE, percentiles = as.numeric(percentiles)),
      AL = data.frame(stringsAsFactors = FALSE, achievementLevel = gsub("Below", "below", as.character(statisticsLevelOut))),
      data.frame()
    )
    if (returnSimpleDoF) {
      # add dofA and dofB if requested
      if (nrow(vec) > 0) {
        vec <- cbind(
          vec,
          data.frame(stringsAsFactors = FALSE, estimateA = resA, estimateAse = SEs[[1]], dofA = JRdfA)
        )
      } else {
        vec <- data.frame(stringsAsFactors = FALSE, estimateA = resA, estimateAse = SEs[[1]], dofA = JRdfA)
      }
    } else {
      # otherwise, no dofA or dofB
      if (nrow(vec) > 0) {
        vec <- cbind(
          vec,
          data.frame(stringsAsFactors = FALSE, estimateA = resA, estimateAse = SEs[[1]])
        )
      } else {
        vec <- data.frame(stringsAsFactors = FALSE, estimateA = resA, estimateAse = SEs[[1]])
      }
    }
    row.names(vec) <- NULL
    if (returnSimpleDoF) {
      percentage <- data.frame(
        stringsAsFactors = FALSE,
        pctA = pctA,
        pctAse = seA,
        dofA = pctDoFA
      )
    } else {
      percentage <- data.frame(
        stringsAsFactors = FALSE,
        pctA = pctA,
        pctAse = seA
      )
    }

    lst <- list(
      results = vec,
      labels = list(
        A = substitute(groupA),
        n0A = nrow2.edsurvey.data.frame(dataA),
        n0B = NA,
        nUsedA = nrow(dA),
        nUsedB = NA
      ),
      percentage = percentage
    )
    if (returnNumberOfPSU) {
      if (getAttributes(dA, "stratumVar") %in% "JK1") {
        lst$labels$nPSUA <- nrow(dA)
      } else {
        lst$labels$nPSUA <- nrow(dA[ , c(getAttributes(dA, "stratumVar"), getAttributes(dA, "psuVar"))])
      }
      lst$labels$nPSUB <- NA
    }
  }

  if (returnVarEstInputs) {
    if ("Level" %in% colnames(varEstInputs$JK)) { # equivalent to type = 'pct' or 'AL'
      if (type == "pct") {
        varEstInputs$JK$Level <- as.numeric(gsub("^P", "", varEstInputs$JK$Level))
        varEstInputs$PV$Level <- as.numeric(gsub("^P", "", varEstInputs$PV$Level))
      }
      if (length(percentiles) <= 1 & length(achievementLevel) <= 1) {
        varEstInputs$JK$Level <- NULL
        varEstInputs$PV$Level <- NULL
      }
    }
    # assign values close to 0 to 0
    varEstInputs$JK$value[which(abs(varEstInputs$JK$value) < (sqrt(.Machine$double.eps) * sqrt(nrow(varEstInputs$JK))))] <- 0
    pctVarEstInputs$JK$value[which(abs(pctVarEstInputs$JK$value) < (sqrt(.Machine$double.eps) * sqrt(nrow(pctVarEstInputs$JK))))] <- 0

    lst <- c(lst, list(
      varEstInputs = varEstInputs,
      pctVarEstInputs = pctVarEstInputs
    ))
  }
  class(lst) <- "gap"
  return(lst)
}

# @title Linking error for PISA 2015 and subsequent DBA v prior PBA
# @param years a vector of two years involved, including 2000, 2003, 2006, 2009, 2012, 2015 or 2018
# @param subject, one of "read", "math", or "scei"
# The linking error for PISA is entirely cross year and one value
# @author Paul Bailey
calLinkingErrorPISA <- function(subject = c("read", "math", "scie"),
                                years = c(2000, 2003, 2006, 2009, 2012, 2015, 2018)) {
  years <- as.numeric(years)
  linkingErrorsFinance <- data.frame(
    subject = rep("fin", 3),
    yearSmall = c(2012, 2015, 2012),
    yearBig = c(2018, 2018, 2015),
    error = c(5.55, 9.37, 5.3309),
    stringsAsFactors = FALSE
  )
  linkingErrorsRead <- data.frame(
    subject = rep("read", 11),
    yearSmall = c(2000, 2003, 2006, 2009, 2012, 2000, 2003, 2006, 2009, 2012, 2015),
    yearBig = c(2015, 2015, 2015, 2015, 2015, 2018, 2018, 2018, 2018, 2018, 2018),
    error = c(6.8044, 5.3907, 6.6064, 3.4301, 5.2535, 4.04, 7.77, 5.24, 3.52, 3.74, 3.93),
    stringsAsFactors = FALSE
  )
  linkingErrorsMath <- data.frame(
    subject = rep("math", 9),
    yearSmall = c(2003, 2006, 2009, 2012, 2003, 2006, 2009, 2012, 2015),
    yearBig = c(2015, 2015, 2015, 2015, 2018, 2018, 2018, 2018, 2018),
    error = c(5.6080, 3.511, 3.7853, 3.5462, 2.80, 3.18, 3.54, 3.34, 2.33),
    stringsAsFactors = FALSE
  )
  linkingErrorsScei <- data.frame(
    subject = rep("scei", 7),
    yearSmall = c(2006, 2009, 2012, 2006, 2009, 2012, 2015),
    yearBig = c(2015, 2015, 2015, 2018, 2018, 2018, 2018),
    error = c(4.4821, 4.5016, 3.9228, 3.47, 3.59, 4.01, 1.51),
    stringsAsFactors = FALSE
  )
  linkingErrorsPISA <- rbind(linkingErrorsFinance, linkingErrorsRead, linkingErrorsMath, linkingErrorsScei)
  sbj <- match.arg(subject)
  years <- years[years %in% c(2000, 2003, 2006, 2009, 2012, 2015, 2018)]
  if (length(years) != 2) {
    stop("Argument ", dQuote("years"), " must have 2 valid years so there are two years to contrast.")
  }
  lep <- linkingErrorsPISA[linkingErrorsPISA$subject == sbj, ]

  res <- expand.grid(yearSmall = years, yearBig = years)
  res <- res[res$yearSmall < res$yearBig, ]
  res <- merge(lep, res, by = c("yearSmall", "yearBig"), all = FALSE)
  if (nrow(res) == 0) {
    return(0)
  }
  if (nrow(res) > 1) {
    stop("Multiple returns.")
  }
  return(res$error)
}

# @title Linking error for NAEP DBA/PBA gap analysis.
# @param X a numeric value estimated from digitally based assessments
# @param subjectI a character value set to "Mathematics" or "Reading" for subject information
# @param gradeI a character value set to "Grade 4" or "Grade 8" as grade level
# @param gapI  a logical indicating if the linking error is a gap
# @param dependentVarI a character value set to one of NAEP math or reading subscales
# @param statElementI a character vector selected from "Mean", "Row Percent", "10", "25", "50", "75", "90",
# "At or above basic", "At or above proficient", "Below basic", "At basic", "At proficient", or "At advanced" indicating the relevant statistic
# @author Huade Huo
calLinkingError <- function(X,
                            subjectI = NULL,
                            gradeI = NULL,
                            dependentVarI = NULL,
                            statElementI = NULL,
                            gapI = FALSE) {
  linkingErrorsRDS <- readRDS(system.file("linkingErrors", "linkingErrors.rds", package = "EdSurvey"))
  linkingErrorRDSi <- linkingErrorsRDS[linkingErrorsRDS$subject %in% subjectI &
    linkingErrorsRDS$dependentVar %in% dependentVarI &
    linkingErrorsRDS$grade %in% gradeI &
    linkingErrorsRDS$statElement %in% statElementI &
    linkingErrorsRDS$gap %in% gapI, ]
  linkingErrorRDSi$statElement <- factor(linkingErrorRDSi$statElement, levels = statElementI)
  linkingErrorRDSi <- linkingErrorRDSi[order(linkingErrorRDSi$statElement), ]
  if (nrow(linkingErrorRDSi) == 0) {
    stop("Linking error is not avaialble.")
  }

  linkingErrori <- linkingErrorRDSi$Aterm * X^2 + linkingErrorRDSi$Bterm * X + linkingErrorRDSi$Cterm
  return(linkingErrori)
}

#' @rdname printGap
#' @title Gap Analysis Printing
#'
#' @description Prints labels and a results vector of a gap analysis.
#' @param x an \code{R} object representing a \code{gap} of class \code{gap} or \code{gapList}
#' @param printPercentage a logical value set to \code{TRUE} to request printing
#'                        of the percentage in the groups. Defaults to \code{TRUE}.
#' @param ... these arguments are not passed anywhere and are included only for compatibility
#' @param use_es_round use the EdSurvey rounding methods for gap
#' @method print gap
#' @author Paul Bailey
#' @aliases print.gapList
#' @export
print.gap <- function(x, ..., printPercentage = TRUE, use_es_round=getOption("EdSurvey_round_output")) {
  if(use_es_round) {
  	x <- es_round(x)
  }
  cat("Call: ")
  print(x$call)
  cat("\nLabels:\n")
  lab <- data.frame(
    stringsAsFactors = FALSE,
    group = c("A", "B"),
    definition = c(deparse(x$labels$A), deparse(x$labels$B)),
    nFullData = c(x$labels$n0A, x$labels$n0B),
    nUsed = c(x$labels$nUsedA, x$labels$nUsedB)
  )
  if (!is.null(x$labels$nPSUA)) {
    lab$nPSU <- c(x$labels$nPSUA, x$labels$nPSUB)
  }
  print(lab, row.names = FALSE)
  cl_print <- match.call()
  if (printPercentage) {
    cat("\nPercentage:\n")
    pct <- data.frame(x$percentage, stringsAsFactors = FALSE)
    res <- prepGapDf(pct)
    res$df <- roundLikePrintCoefmat(res$df, cs=res$cs, tst=res$tst, cl_print)
    print(res$df, row.names = FALSE, cs.ind = res$cs, tst.ind = res$tst, ...)
  }
  cat("\nResults:\n")
  attributes(x$results)$row.names <- rep("", nrow(x$results))
  inds = 1:ncol(x$results)
  tst <- inds[grepl("pValue", colnames(x$results))]
  cs <- inds[!grepl("pValue", colnames(x$results))]
  ind0 <- -1 + which.max(colnames(x$results) == "estimateA")
  tst <- tst[!tst %in% (1:ind0)]
  cs <- cs[!cs %in% (1:ind0)]
  cl_print <- match.call()
  x$results <- roundLikePrintCoefmat(x$results, cs=cs, tst=tst, cl_print)
  print(x$results)
}

roundLikePrintCoefmat <- function(x, cs, tst, cl_print) {
  # round like printCoefmat
  absvec <- abs(x[ , cs])
  absvec <- absvec[absvec > 0 & !is.na(absvec)]
  if("digits" %in% names(cl_print)) {
    digits <- cl_print["digits"]
  } else {
  	digits <- max(3, getOption("digits") - 2)
  }
  if("eps.Pvalue" %in% names(cl_print)) {
    eps.Pvalue <- cl_print["eps.Pvalue"]
  } else {
    eps.Pvalue <- .Machine$double.eps
  }
  if("dig.tst" %in% names(cl_print)) {
  	dig.tst <- cl_print["dig.tst"]
  } else {
  	dig.tst <- max(1, min(5, digits-1))
  }
  new_digits <- max(1, digits - 1+floor(log10(max(absvec))))
  x[,cs] <- format(round(x[ , cs], new_digits), digits=digits)
  x[,tst] <- format.pval(x[ , tst], digits= dig.tst, eps= eps.Pvalue)
  return(x)
}

#' @rdname printGap
#' @method print gapList
#' @export
print.gapList <- function(x, ..., printPercentage = TRUE) {
  cat("gapList\n")
  cat(paste0("Call: "))
  print(x$call)
  cat("\nLabels:\n")
  lab <- data.frame(
    stringsAsFactors = FALSE,
    group = c("A", "B"),
    definition = c(deparse(x$labels$A), deparse(x$labels$B))
  )
  print(lab, row.names = FALSE)
  cl_print <- match.call()
  if (printPercentage) {
    cat("\nPercentage:\n")
    pct <- data.frame(x$percentage, stringsAsFactors = FALSE)
    res <- prepGapDf(pct)
    res$df <- roundLikePrintCoefmat(res$df, cs=res$cs, tst=res$tst, cl_print)
    print(res$df, row.names = FALSE, cs.ind = res$cs, tst.ind = res$tst, ...)
  }
  cat("\nResults:\n")
  res <- prepGapDf(x$results)
  res$df <- roundLikePrintCoefmat(res$df, cs=res$cs, tst=res$tst, cl_print)
  print(res$df)
}

prepGapDf <- function(df) {
  attributes(df)$row.names <- rep("", nrow(df))
  inds = 1:ncol(df)
  tst <- inds[grepl("pValue", colnames(df))]
  cs <- inds[!grepl("pValue", colnames(df))]
  ind0 <- -1 + which.max(colnames(df) %in% c("estimateA", "pctA"))
  tst <- tst[!tst %in% (1:ind0)]
  cs <- cs[!cs %in% (1:ind0)]
  return(list(df=df, tst=tst, cs=cs))
}


# helper that identifies variables in a call
# ccall is the call to parse
# x is an edsurvey.data.frame or light.edsurvey.data.frame
# @author Paul Bailey
parseVars <- function(ccall, x) {
  vars <- c()
  # for each element
  if (typeof(ccall) == "symbol") {
    return(vars)
  }
  for (i in seq_along(ccall)) {
    # if it is a name
    if (inherits(ccall[[i]], "name")) {
      ccall_c <- as.character(ccall[[i]])
      # if it is not in the data and is in the parent.frame, it might be a variable.
      if (ccall_c %in% colnames(x)) { # orignally:: if(ccall_c %in% names(x$data) ) { TOM FINK::dataList update
        if (ccall[[i]] == "%in%" || is.function(ccall[[i]])) {
          # do nothing
        } else {
          vars <- c(vars, ccall_c)
        }
      } # End of if statment: if (ccall_c %in% names(x$data))
    } # end if(class(ccall[[i]]) %in% c("name")) {
    if (inherits(ccall[[i]], "name")) {
      # if this is a call, recursively parse that
      vars <- c(vars, parseVars(ccall[[i]], x))
    } # end of if statment: if class(ccall[[i]]) %in% "call"
  } # end of for loop: i in seq_along(ccall)
  vars
} # End of fucntion: parseVars

# helpter function that calculates "pooled" dof, used for two distinct samples
# @author Trang Nguyen and Paul Bailey
dofCompute <- function(seA, seB, dofA, dofB) {
  return(mapply(function(seA, seB, dofA, dofB) {
    ifelse(seA + seB > 0, (seA^2 + seB^2)^2 / (seA^4 / dofA + seB^4 / dofB), 0)
  }, seA, seB, dofA, dofB))
}

# returns NA when df is 0, but passes Info to pt to avoid a warning
# @author Paul Bailey
pt2 <- function(q, df) {
  df2 <- ifelse(df == 0, Inf, df)
  return(ifelse(df == 0, NA, pt(q, df2)))
}

# match strings to attributes in a case insensitive way
# throwing an error if the string matches more than one achievement level
# or matches none
gmatchAttr <- function(strings, achievementLevelNames) {
  s0 <- strings
  alN0 <- achievementLevelNames
  # lower case everyting to make the function case insensitive
  strings <- tolower(trimws(strings))
  alN <- tolower(achievementLevelNames)
  # find and store the prefix, then remove it from strings
  prefix <- ifelse(grepl("^at or above", strings), substr(s0, 0, nchar("at or above") + 1), "")
  prefix <- ifelse(nchar(prefix) == 0 & grepl("^at ", strings), substr(s0, 0, 3), prefix)
  prefix <- ifelse(grepl("^below", strings), substr(s0, 0, nchar("below") + 1), prefix)
  strings <- trimws(substr(strings, nchar(prefix) + 1, nchar(strings)))
  out <- rep("", length(strings))
  error <- rep("", length(strings))
  for (i in seq_along(out)) {
    for (j in seq_along(alN)) {
      update <- grepl(strings[i], alN[j], fixed = TRUE)
      if (update) {
        if (nchar(out[i]) == 0) {
          # base case, store the output
          out[i] <- ifelse(update, alN0[j], out[i])
        } else {
          # append this level to previously found levels in an error
          # we can then replay all the levels to the user at the end.
          # first...
          # if this is the second match, store the previous
          # one as an error
          if (nchar(error[i]) == 0) {
            error[i] <- sQuote(out[i])
          }
          error[i] <- paste0(c(error[i], sQuote(alN0[j])), collapse = ", ")
          out[i] <- "Error"
        }
      }
    }
  }
  # finalize formatting for errors
  error <- ifelse(nchar(error) > 0, paste0(" matched to ", error), "")
  error <- ifelse(nchar(out) == 0, paste0(" is not one of the achievement Level for the current survey. Choose one from ", pasteItems(achievementLevelNames)), error)
  # if there is an error, format it and throw it
  if (any(nchar(error) > nchar(out))) {
    evars <- nchar(error) > 0
    ferror <- paste0(sQuote(s0[evars]), error[evars], collapse = "; ")
    stop(paste0(ferror, "."))
  }
  out <- trimws(paste(trimws(prefix), trimws(out), sep = " "))
  return(out)
}

getALNames <- function(data, variable) {
  al <- getAttributes(data, "achievementLevels")
  if (is.list(al)) {
    if (variable %in% names(al)) {
      al <- al[[variable]]
    } else {
      # assumes all achievement levels are the same
      al <- al[[1]]
    }
  }
  names(al)
}

addALPrefix <- function(al, als, discrete) {
  if (discrete) {
    if (grepl("below", al, ignore.case = TRUE)) {
      lal <- al
    } else {
      if (!grepl("^at", al, ignore.case = TRUE)) {
        lal <- paste0("At ", al)
      } else {
        lal <- al
      }
    }
  } else {
    if (grepl("^below", al, ignore.case = TRUE)) { # cumulative should not have Below achievementLevel
      return(NULL)
    } else {
      # if it doesn't have the at prefix, add it
      if (!grepl("^at", al, ignore.case = TRUE)) {
        lal <- paste0(
          ifelse(tolower(al) == tolower(als[length(als)]),
            "At ",
            "At or Above "
          ),
          al
        )
      } else {
        lal <- al
      }
    }
  }
  return(lal)
}

Try the EdSurvey package in your browser

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

EdSurvey documentation built on June 27, 2024, 5:10 p.m.