|
| 1 | +#' Build a Complex-Survey Design with Validation and Safe Subpopulations |
| 2 | +#' |
| 3 | +#' @description |
| 4 | +#' A student-friendly wrapper around \code{survey::svydesign()} that takes plain |
| 5 | +#' column names (rather than formulas), validates the inputs with clear error |
| 6 | +#' messages, and -- crucially -- supports subpopulation ("domain") analysis the |
| 7 | +#' \emph{correct} way, by subsetting the design instead of pre-filtering the rows. |
| 8 | +#' |
| 9 | +#' @details |
| 10 | +#' The most common and most consequential mistake in complex-survey analysis is to |
| 11 | +#' \code{filter()} the data to a subgroup \emph{before} building the design. That |
| 12 | +#' discards the design information needed to compute correct standard errors for |
| 13 | +#' the subgroup. The right approach is to build the design on the \strong{full} |
| 14 | +#' sample and then \code{subset()} it. Passing a `subpop` condition here does |
| 15 | +#' exactly that. |
| 16 | +#' |
| 17 | +#' The function also flags strata that contain a single primary sampling unit |
| 18 | +#' (PSU). Such "lonely PSUs" make variance estimation fail unless you set |
| 19 | +#' \code{options(survey.lonely.psu = "adjust")} (or "average"); the function only |
| 20 | +#' \emph{advises} this -- it does not change global options on your behalf. |
| 21 | +#' |
| 22 | +#' @param data A \code{data.frame} containing the survey variables. |
| 23 | +#' @param ids A string naming the cluster / primary sampling unit (PSU) variable |
| 24 | +#' (e.g. \code{"psu"}). Use \code{"0"} or \code{NA} for no clustering. |
| 25 | +#' @param weights A string naming the sampling-weight variable (e.g. |
| 26 | +#' \code{"survey_weight"}). |
| 27 | +#' @param strata A string naming the stratification variable, or \code{NULL} |
| 28 | +#' (the default) for no strata. |
| 29 | +#' @param nest Logical; passed to \code{survey::svydesign()}. If \code{TRUE} |
| 30 | +#' (default), PSU ids are treated as nested within strata. |
| 31 | +#' @param subpop Optional character string giving a logical condition that |
| 32 | +#' defines a subpopulation (e.g. \code{"age >= 20 & sex == 'Female'"}). When |
| 33 | +#' supplied, the design is built on the full sample and then restricted with |
| 34 | +#' \code{survey::subset()}, preserving correct standard errors. |
| 35 | +#' @param verbose Logical; if \code{TRUE} (default), prints the analytic sample |
| 36 | +#' size and any lonely-PSU advisory as messages. |
| 37 | +#' |
| 38 | +#' @return A survey design object of class \code{survey.design2} (the same class |
| 39 | +#' returned by \code{survey::svydesign()}), restricted to the subpopulation |
| 40 | +#' when `subpop` is supplied. Pass it to \code{svytable1()}, \code{svyglm()}, |
| 41 | +#' and the other functions in this package. |
| 42 | +#' |
| 43 | +#' @importFrom survey svydesign |
| 44 | +#' @importFrom stats as.formula |
| 45 | +#' |
| 46 | +#' @export |
| 47 | +#' |
| 48 | +#' @examples |
| 49 | +#' data(nhanes_mortality, package = "svyTable1") |
| 50 | +#' |
| 51 | +#' # Build the design from plain column names. |
| 52 | +#' design <- svydesign_build( |
| 53 | +#' data = nhanes_mortality, |
| 54 | +#' ids = "psu", |
| 55 | +#' strata = "strata", |
| 56 | +#' weights = "survey_weight" |
| 57 | +#' ) |
| 58 | +#' |
| 59 | +#' # Correct subpopulation analysis: build on the full sample, then subset. |
| 60 | +#' design_women <- svydesign_build( |
| 61 | +#' data = nhanes_mortality, |
| 62 | +#' ids = "psu", |
| 63 | +#' strata = "strata", |
| 64 | +#' weights = "survey_weight", |
| 65 | +#' subpop = "sex == 'Female'" |
| 66 | +#' ) |
| 67 | +svydesign_build <- function(data, ids, weights, strata = NULL, |
| 68 | + nest = TRUE, subpop = NULL, verbose = TRUE) { |
| 69 | + |
| 70 | + # --- 1. Input validation --- |
| 71 | + if (!is.data.frame(data)) { |
| 72 | + stop("`data` must be a data.frame.", call. = FALSE) |
| 73 | + } |
| 74 | + if (missing(ids) || !is.character(ids) || length(ids) != 1) { |
| 75 | + stop("`ids` must be a single column name (a string), e.g. \"psu\".", |
| 76 | + call. = FALSE) |
| 77 | + } |
| 78 | + if (missing(weights) || !is.character(weights) || length(weights) != 1) { |
| 79 | + stop("`weights` must be a single column name (a string), e.g. ", |
| 80 | + "\"survey_weight\".", call. = FALSE) |
| 81 | + } |
| 82 | + |
| 83 | + needed <- c(ids, weights) |
| 84 | + if (!is.null(strata)) { |
| 85 | + if (!is.character(strata) || length(strata) != 1) { |
| 86 | + stop("`strata` must be a single column name (a string) or NULL.", |
| 87 | + call. = FALSE) |
| 88 | + } |
| 89 | + needed <- c(needed, strata) |
| 90 | + } |
| 91 | + # ids = "0" is the survey convention for "no clusters"; do not require a column. |
| 92 | + needed <- setdiff(needed, "0") |
| 93 | + missing_cols <- setdiff(needed, names(data)) |
| 94 | + if (length(missing_cols) > 0) { |
| 95 | + stop("Column(s) not found in `data`: ", |
| 96 | + paste(missing_cols, collapse = ", "), ".", call. = FALSE) |
| 97 | + } |
| 98 | + |
| 99 | + w <- data[[weights]] |
| 100 | + if (!is.numeric(w)) { |
| 101 | + stop("The weights column '", weights, "' must be numeric.", call. = FALSE) |
| 102 | + } |
| 103 | + if (anyNA(w)) { |
| 104 | + stop("The weights column '", weights, "' contains missing values. ", |
| 105 | + "Survey weights must be complete.", call. = FALSE) |
| 106 | + } |
| 107 | + if (any(w < 0)) { |
| 108 | + stop("The weights column '", weights, "' contains negative values. ", |
| 109 | + "Survey weights must be non-negative.", call. = FALSE) |
| 110 | + } |
| 111 | + |
| 112 | + # --- 2. Build the design on the FULL sample --- |
| 113 | + ids_f <- stats::as.formula(paste0("~", ids)) |
| 114 | + w_f <- stats::as.formula(paste0("~", weights)) |
| 115 | + strata_f <- if (!is.null(strata)) stats::as.formula(paste0("~", strata)) else NULL |
| 116 | + |
| 117 | + design <- survey::svydesign(ids = ids_f, strata = strata_f, weights = w_f, |
| 118 | + data = data, nest = nest) |
| 119 | + n_full <- nrow(data) |
| 120 | + |
| 121 | + # --- 3. Optional subpopulation via subset (correct for SEs) --- |
| 122 | + if (!is.null(subpop)) { |
| 123 | + if (!is.character(subpop) || length(subpop) != 1) { |
| 124 | + stop("`subpop` must be a single character string giving a logical ", |
| 125 | + "condition, e.g. \"age >= 20\".", call. = FALSE) |
| 126 | + } |
| 127 | + keep <- tryCatch( |
| 128 | + eval(parse(text = subpop), envir = data, enclos = parent.frame()), |
| 129 | + error = function(e) { |
| 130 | + stop("Could not evaluate `subpop` (\"", subpop, "\"): ", |
| 131 | + conditionMessage(e), call. = FALSE) |
| 132 | + } |
| 133 | + ) |
| 134 | + if (!is.logical(keep) || length(keep) != n_full) { |
| 135 | + stop("`subpop` must evaluate to a logical condition over the rows of ", |
| 136 | + "`data` (got length ", length(keep), ").", call. = FALSE) |
| 137 | + } |
| 138 | + keep[is.na(keep)] <- FALSE |
| 139 | + design <- design[which(keep), ] |
| 140 | + if (verbose) { |
| 141 | + message(sprintf( |
| 142 | + "Subpopulation: %d of %d rows retained via subset() (design structure preserved for correct SEs).", |
| 143 | + sum(keep), n_full)) |
| 144 | + } |
| 145 | + } else if (verbose) { |
| 146 | + message(sprintf("Survey design built on %d rows.", n_full)) |
| 147 | + } |
| 148 | + |
| 149 | + # --- 4. Lonely-PSU advisory (does not change global options) --- |
| 150 | + if (verbose && !is.null(strata) && !identical(ids, "0")) { |
| 151 | + psu_per_stratum <- tapply(data[[ids]], data[[strata]], |
| 152 | + function(x) length(unique(x))) |
| 153 | + n_lonely <- sum(psu_per_stratum < 2, na.rm = TRUE) |
| 154 | + if (n_lonely > 0) { |
| 155 | + message(sprintf( |
| 156 | + "Note: %d stratum/strata contain a single PSU. Variance estimation may fail; consider options(survey.lonely.psu = \"adjust\").", |
| 157 | + n_lonely)) |
| 158 | + } |
| 159 | + } |
| 160 | + |
| 161 | + design |
| 162 | +} |
0 commit comments