complexr: Analysis of Complex Survey Data

Overview

complexr is an R package that provides a tidy-oriented framework for the analysis of complex survey data. It supports:

Reading microdata from multiple file formats (CSV, XLSX, SPSS, Stata, RDS)
Constructing and diagnosing stratified and clustered sampling designs
Estimating population means, totals, proportions, ratios, and quantiles with correct variance estimation under complex designs
Deriving new variables from existing survey data
Launching an interactive Shiny application for point-and-click analysis

All estimation functions account for stratification, clustering, and unequal sampling weights following the linearization approach described in Lumley (2010) and Särndal et al. (1992).

1 Conceptual Framework and Notation

This section establishes the statistical notation used throughout the vignette, following the conventions of Gutiérrez et al. (2025).

1.1 Population and Sample

Let $U = \{1, 2, \ldots, N\}$ be the finite population of size $N$, and $s \subset U$ the sample selected under a probability design $p(s)$. For each unit $k \in U$, $y_k$ denotes the value of the variable of interest. The population total and population mean are defined respectively as:

\[ Y = \sum_{k \in U} y_k, \qquad \bar{Y} = \frac{Y}{N}. \]

1.2 Design Weights and Adjusted Weights

The inclusion probability of unit $k$ is $\pi_k = \Pr(k \in s) > 0$. The basic design weight is $d_k = 1/\pi_k$. In practice these weights are modified to incorporate non-response adjustments or calibration to known population totals, yielding adjusted weights $w_k$. Throughout this documentation, $w_k$ refers to the final weights available in the microdata file (variable weight).

1.3 Horvitz–Thompson Estimator

The Horvitz–Thompson (HT) estimator of the population total is (Horvitz and Thompson 1952):

\[ \hat{Y}_{HT} = \sum_{k \in s} d_k\, y_k, \]

and the estimated population size is:

\[ \hat{N}_{HT} = \sum_{k \in s} d_k. \]

When working with adjusted weights $w_k$, the weighted HT estimator takes the form $\hat{Y}_w = \sum_{k \in s} w_k\, y_k$.

1.4 Design-Based Variance

The variance of the HT estimator is estimated as (Särndal et al. 1992):

\[ \hat{V}_p\!\left(\hat{Y}_{HT}\right) = \sum_{k \in s}\sum_{l \in s} \bigl(d_k d_l - d_{kl}\bigr)\, y_k\, y_l, \]

where $d_{kl} = 1/\pi_{kl}$ and $\pi_{kl} = \Pr(k, l \in s)$ are the second-order inclusion probabilities. In practice, equivalent methods such as Taylor linearization or replication (jackknife, bootstrap) are used because they do not require explicit computation of $\pi_{kl}$.

1.5 Stratified Multistage Design

For a design with $H$ strata, $\alpha_h$ primary sampling units (PSUs) in stratum $h$, and $n_{h\alpha}$ observations in PSU $\alpha$, the total estimator is:

\[ \hat{Y}_{HT} = \sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}\, y_{h\alpha k}, \]

where $\omega_{h\alpha k}$ is the adjusted weight of individual $k$ in PSU $\alpha$ of stratum $h$.

1.6 Design Effect (DEFF)

Following Kish (1965), the design effect is defined as the ratio of the estimator variance under the complex design to the variance of the same estimator under simple random sampling (SRS) of the same size:

\[ \widehat{\text{DEFF}} = \frac{\hat{V}_p(\hat{\theta})}{\hat{V}_{\text{SRS}}(\hat{\theta})}. \]

A value $\widehat{\text{DEFF}} > 1$ indicates that the complex design inflates variance relative to SRS, while $\widehat{\text{DEFF}} < 1$ signals an efficiency gain, typically arising from effective stratification.

2 Installation

# Install from GitHub
# install.packages("remotes")
remotes::install_github("stalynGuerrero/complexr")

library(complexr)

3 Simulated Data

The package ships with generate_example_data(), which generates a three-level hierarchical dataset (PSUs → households → individuals) representative of a stratified multistage survey design.

data <- generate_example_data(n_upm = 100, seed = 123)
dplyr::glimpse(data)
#> Rows: 3,438
#> Columns: 15
#> $ strata     <chr> "S3", "S3", "S3", "S3", "S3", "S3", "S3", "S3", "S3", "S3",…
#> $ upm        <chr> "UPM1", "UPM1", "UPM1", "UPM1", "UPM1", "UPM1", "UPM1", "UP…
#> $ hogar_id   <chr> "UPM1_H1", "UPM1_H1", "UPM1_H1", "UPM1_H1", "UPM1_H1", "UPM…
#> $ persona_id <chr> "UPM1_H1_P1", "UPM1_H1_P2", "UPM1_H1_P3", "UPM1_H1_P4", "UP…
#> $ weight     <dbl> 1326.3970, 1326.3970, 1326.3970, 1326.3970, 1326.3970, 851.…
#> $ region     <chr> "South", "South", "South", "South", "South", "Center", "Nor…
#> $ sexo       <chr> "Male", "Female", "Male", "Female", "Female", "Male", "Fema…
#> $ area       <chr> "Urban", "Urban", "Urban", "Urban", "Urban", "Rural", "Rura…
#> $ edad       <dbl> 54, 42, 22, 35, 46, 41, 49, 63, 42, 8, 36, 0, 31, 21, 18, 2…
#> $ educacion  <fct> Higher, Secondary, Higher, Higher, Higher, Higher, Secondar…
#> $ empleo     <fct> Unemployed, Formal, Formal, Informal, Formal, Informal, Inf…
#> $ ingreso_pc <dbl> 1328.0991, 1328.0991, 1328.0991, 1328.0991, 1328.0991, 817.…
#> $ gasto_pc   <dbl> 1100.2257, 1100.2257, 1100.2257, 1100.2257, 1100.2257, 623.…
#> $ pobre      <int> 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,…
#> $ ingreso2   <dbl> 1143.2572, NA, 1363.3893, 1477.1621, 1583.9139, NA, 663.904…

3.1 Variables

Variable	Type	Description
`strata`	character	Stratum identifier ($h = 1,\ldots,H$)
`upm`	character	Primary sampling unit $\alpha$ within stratum $h$
`hogar_id`	character	Household identifier
`persona_id`	character	Individual identifier $k$
`weight`	numeric	Adjusted weight $w_k$ (inverse inclusion probability, calibrated)
`region`	character	Estimation domain: North / Center / South
`sexo`	character	Sex: Male / Female
`area`	character	Area: Urban / Rural
`edad`	numeric	Age in years
`educacion`	factor	Education level: Primary / Secondary / Higher
`empleo`	factor	Employment status: Formal / Informal / Unemployed
`ingreso_pc`	numeric	Per-capita household income ($y_k$)
`gasto_pc`	numeric	Per-capita household expenditure
`pobre`	numeric	Binary poverty indicator: $y_k \in \{0, 1\}$
`ingreso2`	numeric	Auxiliary income variable (10 % missing values)

3.2 Hierarchy and Internal Consistency

The simulation enforces the following constraints:

Weights are constant within households: $w_{h\alpha k} = w_{h\alpha}$ for all $k$ in the household.
Per-capita income is constant within households: $y_{h\alpha k} = Y_{h\alpha} / N_{h\alpha}$.
Individuals younger than 5 years have NA for education.

Household income follows a hierarchical gamma model with PSU and household random effects:

\[ Y_{h\alpha} \sim \text{Gamma}\!\left(\alpha_0,\; \beta_0 \cdot \exp(u_{h\alpha} + v_{h\alpha k})\right) \]

where $u_{h\alpha} \sim N(0, 0.09)$ is the PSU effect and $v_{h\alpha k} \sim N(0, 0.04)$ is the household effect.

4 Reading External Data

To load your own microdata use read_survey_data(). Supported formats are detected automatically from the file extension.

# CSV
data <- read_survey_data("survey.csv")

# SPSS
data <- read_survey_data("survey.sav")

# Stata
data <- read_survey_data("survey.dta")

# Excel
data <- read_survey_data("survey.xlsx")

The function returns the data as a tibble and attaches metadata attributes: source_path, source_format, n_rows, n_cols.

5 Deriving Variables

mutate_survey_data() creates new variables from a named list of one-sided formulas, evaluated sequentially in the environment of the data frame.

data <- mutate_survey_data(
  data,
  definitions = list(
    log_ingreso = ~ log(ingreso_pc + 1),
    ratio_gasto = ~ gasto_pc / ingreso_pc
  )
)

dplyr::select(data, ingreso_pc, log_ingreso, ratio_gasto) |> head(4)
#> # A tibble: 4 × 3
#>   ingreso_pc log_ingreso ratio_gasto
#>        <dbl>       <dbl>       <dbl>
#> 1      1328.        7.19       0.828
#> 2      1328.        7.19       0.828
#> 3      1328.        7.19       0.828
#> 4      1328.        7.19       0.828

6 Survey Design

6.1 Building a Design Object

as_survey_design_tbl() wraps survey::svydesign() and returns a tbl_svy object compatible with the srvyr / survey ecosystem (Lumley 2010).

design <- as_survey_design_tbl(
  data    = data,
  weight  = "weight",
  strata  = "strata",
  cluster = "upm",
  nest    = TRUE
)

class(design)
#> [1] "tbl_svy"        "survey.design2" "survey.design"

Supported configurations:

Configuration	Arguments
Simple random sampling (SRS)	`weight` only
Stratified	`weight` + `strata`
Clustered (single-stage)	`weight` + `cluster`
Stratified multistage	`weight` + `strata` + `cluster`
With finite population correction	any of the above + `fpc`

The function validates that weights $w_k$ are strictly positive ($w_k > 0$) and free of missing values, and optionally checks that PSUs are not shared across strata (check_psu = TRUE).

Note: When a stratum $h$ contains only one PSU ($\alpha_h = 1$), estimation of $\hat{V}_p$ by Taylor linearization is undefined. The function automatically sets options(survey.lonely.psu = "adjust") to use the conservative centered-at-stratum-mean approximation (Cochran 1977).

6.2 Diagnosing the Design

describe_survey_design(design)
#> # A tibble: 1 × 7
#>   n_obs n_strata n_clusters weight_min weight_max weight_mean weight_cv
#>   <int>    <int>      <int>      <dbl>      <dbl>       <dbl>     <dbl>
#> 1  3438        5        100       205.      1500.        848.     0.444

The diagnostic table reports:

Column	Description
`n_obs`	Total sample size $n = \lvert s \rvert$
`n_strata`	Number of strata $H$
`n_clusters`	Total number of PSUs $\sum_h \alpha_h$
`weight_min`	$\min_{k \in s} w_k$
`weight_max`	$\max_{k \in s} w_k$
`weight_mean`	$\bar{w} = \hat{N}_w / n$
`weight_cv`	$CV(w) = s_w / \bar{w}$

7 Estimation

All estimators are computed with estimate_survey(). The function returns a tibble with the following columns:

Column	Description
`variable`	Name of the target variable
`estimator`	Type of estimator
`estimate`	Point estimate $\hat{\theta}$
`se`	Standard error $ee(\hat{\theta}) = \sqrt{\hat{V}_p(\hat{\theta})}$
`cv`	Coefficient of variation $CV = ee(\hat{\theta})/\hat{\theta}$
`deff`	Design effect $\widehat{\text{DEFF}}$
`lci`	Lower confidence bound
`uci`	Upper confidence bound
`quality`	Precision label based on $CV$

Precision labels (based on the coefficient of variation):

$CV$	Label
$< 5\%$	Very high precision
$5\%$–$10\%$	High precision
$10\%$–$20\%$	Acceptable precision
$20\%$–$30\%$	Use with caution
$\geq 30\%$	Low precision

7.1 Population Mean

The weighted total and estimated population size are:

\[ \hat{Y}_w = \sum_{k \in s} w_k\, y_k, \qquad \hat{N}_w = \sum_{k \in s} w_k. \]

The weighted mean (Horvitz–Thompson ratio estimator) is:

\[ \bar{y}_w = \frac{\hat{Y}_w}{\hat{N}_w} = \frac{\displaystyle\sum_{k \in s} w_k\, y_k} {\displaystyle\sum_{k \in s} w_k}. \]

Its variance is estimated by Taylor linearization (Särndal et al. 1992):

\[ \hat{V}_p\!\left(\bar{y}_w\right) = \frac{1}{\hat{N}_w^2}\,\hat{V}_p\!\left(\hat{Y}_w\right). \]

r_mean <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean"
)
r_mean
#> # A tibble: 1 × 9
#>   variable   estimator estimate    se     cv  deff   lci   uci quality          
#>   <chr>      <chr>        <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>            
#> 1 ingreso_pc mean         1579.  71.6 0.0453  8.05 1439. 1720. Very high precis…

7.2 Population Total

Under a stratified multistage design, the HT total estimator is:

\[ \hat{Y}_w = \sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h}\sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}\, y_{h\alpha k}. \]

Its variance is estimated stratum by stratum:

\[ \hat{V}_p\!\left(\hat{Y}_w\right) = \sum_{h=1}^{H} \hat{V}_{p,h}\!\left(\hat{Y}_{w,h}\right), \]

where $\hat{V}_{p,h}$ is computed within each stratum $h$ using the deviations of PSU totals from their stratum mean.

r_total <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "total"
)
r_total
#> # A tibble: 1 × 9
#>   variable   estimator    estimate         se     cv  deff    lci    uci quality
#>   <chr>      <chr>           <dbl>      <dbl>  <dbl> <dbl>  <dbl>  <dbl> <chr>  
#> 1 ingreso_pc total     4607014476. 261136568. 0.0567  12.6 4.10e9 5.12e9 High p…

7.3 Proportions

7.3.1 Binary Variable

For an indicator variable $y_k \in \{0, 1\}$, the population proportion $\pi$ is estimated as (Heeringa et al. 2017):

\[ \hat{p} = \frac{\displaystyle\sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}\, I(y_k = 1)} {\displaystyle\sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}} = \frac{\hat{N}_1}{\hat{N}_w}. \]

The variance of $\hat{p}$ is approximated by Taylor linearization:

\[ \hat{V}_p(\hat{p}) \;\dot{=}\; \frac{\hat{V}_p(\hat{N}_1) + \hat{p}^2\,\hat{V}_p(\hat{N}_w) - 2\hat{p}\,\widehat{\text{cov}}(\hat{N}_1, \hat{N}_w)} {\hat{N}_w^2}. \]

r_pobre <- estimate_survey(
  design    = design,
  variable  = "pobre",
  estimator = "prop"
)
r_pobre
#> # A tibble: 1 × 9
#>   variable estimator estimate      se     cv  deff   lci   uci quality          
#>   <chr>    <chr>        <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>            
#> 1 pobre    prop         0.259 0.00936 0.0361  1.57 0.241 0.277 Very high precis…

7.3.2 Multinomial Categorical Variable

For a variable with categories $\mathcal{K} = \{k_1, k_2, \ldots\}$, the proportion for category $k$ is:

\[ \hat{p}_k = \frac{\displaystyle\sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{i=1}^{n_{h\alpha}} \omega_{h\alpha i}\, I(y_i = k)} {\displaystyle\sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{i=1}^{n_{h\alpha}} \omega_{h\alpha i}} = \frac{\hat{N}_k}{\hat{N}_w}. \]

The function automatically constructs the indicator $I(y_i = k)$ for each category:

r_empleo <- estimate_survey(
  design    = design,
  variable  = "empleo",
  estimator = "prop"
)
r_empleo
#> # A tibble: 3 × 10
#>   variable estimator estimate      se     cv  deff empleo      lci   uci quality
#>   <chr>    <chr>        <dbl>   <dbl>  <dbl> <dbl> <chr>     <dbl> <dbl> <chr>  
#> 1 empleo   prop         0.331 0.00815 0.0246  1.03 Formal    0.315 0.347 Very h…
#> 2 empleo   prop         0.341 0.00903 0.0265  1.25 Informal  0.323 0.359 Very h…
#> 3 empleo   prop         0.328 0.00866 0.0264  1.17 Unemploy… 0.311 0.345 Very h…

7.4 Ratio Estimator

The ratio estimator of two population totals is (Cochran 1977):

\[ \hat{R} = \frac{\hat{Y}_w}{\hat{X}_w} = \frac{\displaystyle\sum_{k \in s} w_k\, y_k} {\displaystyle\sum_{k \in s} w_k\, x_k}. \]

Variance is estimated by first-order Taylor linearization:

\[ \hat{V}_p(\hat{R}) \approx \frac{1}{\hat{X}_w^2}\, \hat{V}_p\!\left(\hat{Y}_w - \hat{R}\,\hat{X}_w\right). \]

Numeric / Numeric:

r_ratio <- estimate_survey(
  design      = design,
  estimator   = "ratio",
  numerator   = "ingreso_pc",
  denominator = "gasto_pc"
)
r_ratio
#> # A tibble: 1 × 9
#>   variable          estimator estimate      se      cv  deff   lci   uci quality
#>   <chr>             <chr>        <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <chr>  
#> 1 ingreso_pc_over_… ratio         1.34 0.00579 0.00432    NA  1.33  1.35 Very h…

Categorical / Categorical — ratio of Formal to Informal workers ($\hat{N}_{\text{Formal}} / \hat{N}_{\text{Informal}}$):

r_ratio_cat <- estimate_survey(
  design          = design,
  estimator       = "ratio",
  numerator       = "empleo",
  denominator     = "empleo",
  ratio_num_level = "Formal",
  ratio_den_level = "Informal"
)
r_ratio_cat
#> # A tibble: 1 × 9
#>   variable           estimator estimate     se     cv  deff   lci   uci quality 
#>   <chr>              <chr>        <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>   
#> 1 empleo_over_empleo ratio        0.970 0.0429 0.0442    NA 0.886  1.05 Very hi…

Numeric / Categorical — average income among formal workers, equivalent to $\hat{Y}_{\text{income}} / \hat{N}_{\text{Formal}}$:

r_ratio_mix <- estimate_survey(
  design          = design,
  estimator       = "ratio",
  numerator       = "ingreso_pc",
  denominator     = "empleo",
  ratio_den_level = "Formal"
)
r_ratio_mix
#> # A tibble: 1 × 9
#>   variable             estimator estimate    se     cv  deff   lci   uci quality
#>   <chr>                <chr>        <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>  
#> 1 ingreso_pc_over_emp… ratio        4773.  242. 0.0507    NA 4299. 5247. High p…

7.5 Quantiles

Quantiles are derived from the weighted empirical cumulative distribution function (Woodruff 1952):

\[ \hat{F}_w(t) = \frac{\displaystyle\sum_{k \in s} w_k\, I(y_k \le t)} {\displaystyle\sum_{k \in s} w_k} = \frac{\hat{N}(y \le t)}{\hat{N}_w}. \]

The $p$-th order quantile is defined as:

\[ \hat{q}_p = \inf\bigl\{t : \hat{F}_w(t) \ge p\bigr\}. \]

Confidence intervals are computed using the Woodruff linearization method, which transforms the problem to the scale of the cumulative proportion:

\[ IC_p[\hat{q}_p] = \left\{t : \hat{F}_w(t) \in \left[p \pm t_{1-\alpha/2,\,df}\; ee(\hat{F}_w(t))\right]\right\}. \]

r_quant <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "quantile",
  probs     = c(0.10, 0.25, 0.50, 0.75, 0.90)
)
r_quant
#> # A tibble: 5 × 10
#>   variable   estimator quantile estimate    se     cv  deff   lci   uci quality 
#>   <chr>      <chr>        <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>   
#> 1 ingreso_pc quantile      0.1      499.  28.0 0.0562    NA  444.  554. High pr…
#> 2 ingreso_pc quantile      0.25     736.  32.3 0.0438    NA  673.  800. Very hi…
#> 3 ingreso_pc quantile      0.5     1127.  42.0 0.0373    NA 1045. 1209. Very hi…
#> 4 ingreso_pc quantile      0.75    1895. 106.  0.0561    NA 1686. 2103. High pr…
#> 5 ingreso_pc quantile      0.9     3196. 221.  0.0691    NA 2763. 3629. High pr…

8 Domain Estimation

In household surveys it is common to estimate parameters for subpopulations or domains $U_d \subset U$. The weighted ratio estimator in domain $d$ is:

\[ \bar{y}_{w,d} = \frac{\displaystyle\sum_{k \in s} w_k\, y_k\, I(k \in U_d)} {\displaystyle\sum_{k \in s} w_k\, I(k \in U_d)} = \frac{\hat{Y}_{w,d}}{\hat{N}_{w,d}}. \]

Variance estimation is carried out over the full sample $s$, preserving the design structure and avoiding subsetting bias (Lumley 2010).

8.1 Domain Design Effect

The design effect for domain $d$ is defined analogously to the overall DEFF (Kish 1965):

\[ \widehat{\text{DEFF}}_d = \frac{\hat{V}_p(\hat{\theta}_d)} {\hat{V}_{\text{SRS}}(\hat{\theta}_d)}, \]

where $\hat{V}_{\text{SRS}}(\hat{\theta}_d)$ is the variance that would be obtained under SRS restricted to domain $d$. For the domain mean, this simplifies to:

\[ \hat{V}_{\text{SRS}}(\bar{y}_{w,d}) = \left(1 - \frac{n_d}{N_d}\right)\frac{S_{y,d}^2}{n_d}, \]

with $n_d = \sum_{k \in s} I(k \in U_d)$ the domain sample size, $N_d \approx \hat{N}_{w,d}$ the estimated domain population size, and $S_{y,d}^2$ the unweighted within-domain sample variance. Values $\widehat{\text{DEFF}}_d > 1$ indicate that clustering or unequal weighting inflate variance even within the domain.

The by argument controls estimation domains:

r_region <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean",
  by        = "region"
)
r_region
#> # A tibble: 3 × 10
#>   variable   estimator region estimate    se     cv  deff   lci   uci quality   
#>   <chr>      <chr>     <chr>     <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>     
#> 1 ingreso_pc mean      Center    1597.  88.6 0.0555  4.03 1423. 1770. High prec…
#> 2 ingreso_pc mean      North     1616.  85.2 0.0527  3.62 1449. 1783. High prec…
#> 3 ingreso_pc mean      South     1523.  71.1 0.0467  2.89 1384. 1662. Very high…

Multiple domain variables are supported (crossed domains $U_{d_1} \cap U_{d_2}$):

r_region_area <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean",
  by        = c("region", "area")
)
r_region_area
#> # A tibble: 6 × 11
#> # Groups:   region [3]
#>   variable   estimator region area  estimate    se     cv  deff   lci   uci
#>   <chr>      <chr>     <chr>  <chr>    <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1 ingreso_pc mean      Center Rural    1642. 109.  0.0663  2.91 1429. 1856.
#> 2 ingreso_pc mean      Center Urban    1541. 110.  0.0713  3.44 1325. 1756.
#> 3 ingreso_pc mean      North  Rural    1687. 115.  0.0679  2.76 1462. 1911.
#> 4 ingreso_pc mean      North  Urban    1545.  90.9 0.0588  2.57 1367. 1723.
#> 5 ingreso_pc mean      South  Rural    1512.  91.6 0.0606  2.39 1332. 1691.
#> 6 ingreso_pc mean      South  Urban    1536.  89.0 0.0579  2.29 1362. 1711.
#> # ℹ 1 more variable: quality <chr>

Proportions by domain:

r_pobre_region <- estimate_survey(
  design    = design,
  variable  = "pobre",
  estimator = "prop",
  by        = "region"
)
r_pobre_region
#> # A tibble: 3 × 10
#>   variable estimator region estimate     se     cv  deff   lci   uci quality    
#>   <chr>    <chr>     <chr>     <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>      
#> 1 pobre    prop      Center    0.250 0.0155 0.0621  1.42 0.220 0.281 High preci…
#> 2 pobre    prop      North     0.252 0.0132 0.0523  1.11 0.227 0.278 High preci…
#> 3 pobre    prop      South     0.275 0.0148 0.0537  1.24 0.246 0.304 High preci…

9 Formatting Results

format_results_table() rounds numeric columns, computes missing $CV$ or confidence intervals, and ensures the output always contains columns estimate, se, cv, lci, and uci.

format_results_table(r_region, digits = 3)
#> # A tibble: 3 × 10
#>   variable   estimator region estimate    se    cv  deff   lci   uci quality    
#>   <chr>      <chr>     <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      
#> 1 ingreso_pc mean      Center    1597.  88.6 0.056  4.03 1423. 1770. High preci…
#> 2 ingreso_pc mean      North     1616.  85.2 0.053  3.62 1449. 1783. High preci…
#> 3 ingreso_pc mean      South     1523.  71.1 0.047  2.89 1384. 1662. Very high …

10 Visualisation

plot_results_bar() generates a ggplot2 bar chart with error bars representing the confidence interval $[\hat{\theta} - t\, ee(\hat{\theta}),\; \hat{\theta} + t\, ee(\hat{\theta})]$. Domain variables are detected automatically as any column not in the standard output columns.

plot_results_bar(r_region)

$Per-capita income by region — $\bar{y}_{w,d}$ with 95% CI$

Per-capita income by region — $\bar{y}_{w,d}$ with 95% CI

plot_results_bar(r_region_area)

Per-capita income by region and area, 95% CI

Proportion plots automatically constrain the y-axis to $[0, 1]$:

plot_results_bar(r_pobre_region)

$Poverty rate by region — $\hat{p}_d$ with 95% CI$

Poverty rate by region — $\hat{p}_d$ with 95% CI

11 Shiny Application

The package includes a full interactive Shiny application covering the complete analysis pipeline:

Data — upload CSV / RDS / XLSX or load built-in example data
Survey design — choose SRS, stratified, or multistage cluster; select weight $w_k$, stratum $h$, and PSU $\alpha$ variables; view MathJax theory
Estimation — select estimator ($\bar{y}_w$, $\hat{Y}_w$, $\hat{p}_k$, $\hat{R}$, $\hat{q}_p$), target variable, and domains; execute and inspect results with $CV$-based precision indicators

Launch the application with:

ComplexSurvey_app()

12 Complete Workflow Example

library(complexr)

# 1. Generate / load data
data <- generate_example_data(n_upm = 100, seed = 2024)

# 2. Derive new variables
data <- mutate_survey_data(
  data,
  definitions = list(
    log_ingreso = ~ log(ingreso_pc + 1)
  )
)

# 3. Build survey design (stratified multistage)
design <- as_survey_design_tbl(
  data    = data,
  weight  = "weight",
  strata  = "strata",
  cluster = "upm",
  nest    = TRUE
)

# 4. Diagnose: N_hat, H, PSUs, CV(w)
describe_survey_design(design)
#> # A tibble: 1 × 7
#>   n_obs n_strata n_clusters weight_min weight_max weight_mean weight_cv
#>   <int>    <int>      <int>      <dbl>      <dbl>       <dbl>     <dbl>
#> 1  3371        5        100       208.      1498.        850.     0.439

# 5. Estimate domain mean and format
res <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean",
  by        = c("region", "area")
)

format_results_table(res, digits = 2)
#> # A tibble: 6 × 11
#>   variable estimator region area  estimate    se    cv  deff   lci   uci quality
#>   <chr>    <chr>     <chr>  <chr>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>  
#> 1 ingreso… mean      Center Rural    1559. 103.   0.07  2.13 1357. 1762. High p…
#> 2 ingreso… mean      Center Urban    1482.  85.7  0.06  3.03 1314. 1650. High p…
#> 3 ingreso… mean      North  Rural    1540. 106.   0.07  2.76 1331. 1748. High p…
#> 4 ingreso… mean      North  Urban    1512.  71.0  0.05  1.71 1373. 1652. Very h…
#> 5 ingreso… mean      South  Rural    1520.  86.8  0.06  2.65 1349. 1690. High p…
#> 6 ingreso… mean      South  Urban    1474.  77.3  0.05  1.89 1322. 1625  High p…

References

Cochran, William G. 1977. Sampling Techniques. 3rd ed. John Wiley & Sons.

Gutiérrez, Andrés, Stalyn Guerrero, Cristian Téllez, and Giovany Babativa. 2025. Análisis de Encuestas Con R. https://psirusteam.github.io/2021ASDA/.

Heeringa, Steven G., Brady T. West, and Patricia A. Berglund. 2017. Applied Survey Data Analysis. 2nd ed. Chapman; Hall/CRC. https://doi.org/10.1201/9781315153278.

Horvitz, Daniel G., and Donovan J. Thompson. 1952. “A Generalization of Sampling Without Replacement from a Finite Universe.” Journal of the American Statistical Association 47 (260): 663–85. https://doi.org/10.1080/01621459.1952.10483446.

Kish, Leslie. 1965. Survey Sampling. John Wiley & Sons.

Lumley, Thomas. 2010. Complex Surveys: A Guide to Analysis Using R. John Wiley & Sons. https://doi.org/10.1002/9780470580066.

Särndal, Carl-Erik, Bengt Swensson, and Jan Wretman. 1992. Model Assisted Survey Sampling. Springer-Verlag. https://doi.org/10.1007/978-1-4612-4378-6.

Woodruff, Ralph S. 1952. “Confidence Intervals for Medians and Other Position Measures.” Journal of the American Statistical Association 47 (260): 635–46. https://doi.org/10.1080/01621459.1952.10483443.