---
title: "Basic Workflow with svySE"
author: "Luis Burgos"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Basic Workflow with svySE}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center"
)
```

# Introduction

`svySE` provides a structured workflow for producing indicator tables from
complex survey data.

The package separates two complementary tasks:

1. `svySE_calc()` calculates weighted estimates and sampling errors using
   sampling weights and, when available, strata and cluster variables.
2. `svySE_simple()` calculates unweighted frequencies and percentages from the
   observed sample without using the survey design.

Results generated by either function can be exported to `.xlsx` files with
`svySE_xlsx()`. The export function can process one result or consolidate
several analyses generated from different datasets, indicators, weights, or
function calls.

The package is built on top of `survey` for design-based estimation and uses
`openxlsx` for workbook generation.

# Package workflow

| Step | Function | Purpose |
|------|----------|---------|
| 1 | `svySE_cfg()` | Configure estimation options |
| 2A | `svySE_calc()` | Calculate weighted estimates and sampling errors |
| 2B | `svySE_simple()` | Calculate unweighted frequencies and percentages |
| 3 | `svySE_xlsx()` | Export one or multiple results to `.xlsx` files |

The two calculation functions are intentionally separated. This avoids
performing unnecessary calculations when only sampling errors or only
descriptive tables are required.

# Example data

The following simulated dataset contains:

- one grouping variable;
- one stratification variable;
- one cluster variable;
- one sampling weight;
- one division variable;
- two binary indicators.

```{r example-data}
library(svySE)

set.seed(123)

df <- data.frame(
  dept = rep(c("A", "B", "C"), each = 50),
  strata = rep(c("S1", "S2", "S3"), each = 50),
  cluster = rep(seq_len(30), each = 5),
  service = rep(c("S1", "S2"), length.out = 150),
  weight = runif(150, 10, 50),
  ind_1 = sample(c(0, 1), 150, replace = TRUE),
  ind_2 = sample(c(0, 1), 150, replace = TRUE),
  stringsAsFactors = FALSE
)

head(df)
```

# Configure sampling error estimation

`svySE_cfg()` defines the settings used by `svySE_calc()`.

```{r configuration}
cfg <- svySE_cfg(
  estimator = "prop",
  variance = "taylor",
  lonely_psu = "adjust",
  conf_level = 0.95,
  target = 1,
  valid_values = c(0, 1),
  truncate_lower_ci = TRUE,
  pct_mult = 100,
  deff = TRUE,
  cv = TRUE,
  na_rm = TRUE,
  ci_method = "wald",
  ci_df = NULL
)

cfg
```

The principal options are:

| Argument | Purpose |
|----------|---------|
| `estimator` | Selects the estimator used in the analysis |
| `variance` | Defines the variance estimation method |
| `lonely_psu` | Specifies how strata containing one PSU are handled |
| `conf_level` | Defines the confidence level |
| `target` | Identifies the indicator category of interest |
| `valid_values` | Defines the valid indicator values |
| `truncate_lower_ci` | Truncates negative lower confidence limits at zero |
| `pct_mult` | Converts proportions to the selected percentage scale |
| `deff` | Enables or disables design-effect calculation |
| `cv` | Enables or disables coefficient-of-variation calculation |
| `na_rm` | Controls the treatment of missing values |
| `ci_method` | Selects the confidence interval method for proportions: `"wald"` or `"xlogit"` |
| `ci_df` | Sets the degrees of freedom of the `"xlogit"` interval (`NULL` = design) |

# Calculate sampling errors

`svySE_calc()` calculates weighted estimates and measures of sampling
precision.

```{r sampling-errors}
res_error <- svySE_calc(
  data = df,
  indicators = c("ind_1", "ind_2"),
  group_vars = "dept",
  group_labels = "Department",
  strata = "strata",
  cluster = "cluster",
  weight = "weight",
  division = NULL,
  div_weight = NULL,
  cfg = cfg,
  verbose = FALSE
)

res_error
```

The returned object has class `"svySE_result"`.

```{r sampling-error-class}
class(res_error)
```

Results are organized by indicator and division.

```{r inspect-error-object}
names(res_error$results)
names(res_error$results$ind_1$error)
```

The national and grouped results for the first indicator can be inspected with:

```{r inspect-error-table}
res_error$results$ind_1$error$TOTAL
```

The sampling error table contains:

| Column | Description |
|--------|-------------|
| `est_abs` | Weighted absolute estimate |
| `est_pct` | Weighted percentage estimate |
| `se_abs` | Standard error of the absolute estimate |
| `se_pct` | Standard error of the percentage |
| `ci_l_abs` | Lower confidence limit for the absolute estimate |
| `ci_l_pct` | Lower confidence limit for the percentage |
| `ci_u_abs` | Upper confidence limit for the absolute estimate |
| `ci_u_pct` | Upper confidence limit for the percentage |
| `cv` | Coefficient of variation |
| `deff` | Design effect |
| `n_unw` | Unweighted number of target cases |

# Confidence Intervals for Proportions

Starting with version 0.3.0, `svySE` offers two methods to build the
confidence intervals of proportions (`ci_l_pct` and `ci_u_pct`). The method is
selected with the `ci_method` argument of `svySE_cfg()`.

| Method | `ci_method` | Interval | Critical value |
|--------|-------------|----------|----------------|
| Wald (default) | `"wald"` | `p ± z × SE` | Normal quantile |
| Logit | `"xlogit"` | `expit(logit(p) ± t × SE / (p(1 − p)))` | `t` quantile with the design degrees of freedom |

Both methods use exactly the same estimate and the same Taylor-linearization
standard error. Only the construction of the interval changes: estimates,
standard errors, coefficients of variation, design effects, totals, and counts
are identical under both methods.

The logit interval is computed in four steps:

1. The proportion is transformed to the logit scale: `logit(p) = log(p / (1 − p))`.
2. The standard error is transferred to that scale with the delta method: `SE_logit = SE / (p(1 − p))`.
3. The limits are computed as `logit(p) ± t(1 − α/2, df) × SE_logit`.
4. The limits are transformed back to the 0–1 scale with `expit(x) = 1 / (1 + exp(−x))`.

Main properties of the logit interval:

- limits always remain within `[0, 1]`, without truncation;
- the interval is asymmetric around `p` when the proportion is close to 0 or 1;
- intervals of complementary categories are complementary: the interval of `1 − p` is `[1 − upper, 1 − lower]`;
- the critical value uses a `t` distribution with the design degrees of freedom (number of PSUs minus number of strata; without clusters, number of observations minus number of strata);
- domains use the degrees of freedom of the full design, not those of the domain.

The logit interval is recommended for small or large proportions and for
domains with few cases, where the Wald interval may produce limits outside
`[0, 1]`.

## Use the logit interval

```{r xlogit-calc}
cfg_xlogit <- svySE_cfg(
  estimator = "prop",
  ci_method = "xlogit"
)

res_xlogit <- svySE_calc(
  data = df,
  indicators = c("ind_1", "ind_2"),
  group_vars = "dept",
  group_labels = "Department",
  strata = "strata",
  cluster = "cluster",
  weight = "weight",
  cfg = cfg_xlogit,
  verbose = FALSE
)

res_xlogit$results$ind_1$error$TOTAL
```

## Compare both methods

```{r xlogit-compare}
wald <- res_error$results$ind_1$error$TOTAL
xlogit <- res_xlogit$results$ind_1$error$TOTAL

data.frame(
  dept = wald$dept,
  est_pct = wald$est_pct,
  se_pct = wald$se_pct,
  wald_lower = wald$ci_l_pct,
  wald_upper = wald$ci_u_pct,
  xlogit_lower = xlogit$ci_l_pct,
  xlogit_upper = xlogit$ci_u_pct
)
```

## Degrees of freedom

By default (`ci_df = NULL`), the critical value uses the design degrees of
freedom, which can be inspected with `survey::degf()`:

```{r xlogit-degf}
design <- survey::svydesign(
  ids = ~cluster,
  strata = ~strata,
  weights = ~weight,
  data = df,
  nest = TRUE
)

survey::degf(design)
```

The degrees of freedom can also be set explicitly with `ci_df`. Use
`ci_df = Inf` to apply the normal quantile.

```{r xlogit-df}
cfg_xlogit_df <- svySE_cfg(
  estimator = "prop",
  ci_method = "xlogit",
  ci_df = 30
)

cfg_xlogit_z <- svySE_cfg(
  estimator = "prop",
  ci_method = "xlogit",
  ci_df = Inf
)

cfg_xlogit_df
```

When `division` is used, each division is estimated with its own survey design,
so its degrees of freedom correspond to the records of that division. Use
`ci_df` if a common value is required.

## Complementary categories

With the logit interval, the interval of the complementary category is obtained
directly from the interval of the target category.

```{r xlogit-complement}
res_0 <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  strata = "strata",
  cluster = "cluster",
  weight = "weight",
  cfg = svySE_cfg(estimator = "prop", target = 0, ci_method = "xlogit"),
  verbose = FALSE
)

tab_1 <- res_xlogit$results$ind_1$error$TOTAL
tab_0 <- res_0$results$ind_1$error$TOTAL

all.equal(tab_0$ci_l_pct, 100 - tab_1$ci_u_pct)
all.equal(tab_0$ci_u_pct, 100 - tab_1$ci_l_pct)
```

## Special cases

| Situation | Logit interval |
|-----------|----------------|
| Proportion equal to 0 or 1 | Degenerate interval `[p, p]`, as with the Wald method |
| Standard error equal to 0 | Degenerate interval `[p, p]` |
| Missing estimate or standard error | `NA` limits |
| Degrees of freedom not positive | `NA` limits |
| Groups without observations | `NA` row, as with the Wald method |
| `estimator` other than `"prop"` | Not available; `svySE_cfg()` returns an error |

The confidence intervals of absolute estimates (`ci_l_abs`, `ci_u_abs`) always
use the Wald method.

# Supported survey designs

`svySE_calc()` can be used with different combinations of design variables.

| Design | `strata` | `cluster` |
|--------|----------|-----------|
| Weight only | `NULL` | `NULL` |
| Stratified | Variable name | `NULL` |
| Clustered | `NULL` | Variable name |
| Stratified and clustered | Variable name | Variable name |

## Weight-only design

```{r design-weight-only}
res_weight <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  strata = NULL,
  cluster = NULL,
  weight = "weight",
  cfg = cfg,
  verbose = FALSE
)

res_weight
```

## Stratified design

```{r design-stratified}
res_strata <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  strata = "strata",
  cluster = NULL,
  weight = "weight",
  cfg = cfg,
  verbose = FALSE
)

res_strata
```

## Clustered design

```{r design-clustered}
res_cluster <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  strata = NULL,
  cluster = "cluster",
  weight = "weight",
  cfg = cfg,
  verbose = FALSE
)

res_cluster
```

## Stratified and clustered design

```{r design-complex}
res_complex <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  strata = "strata",
  cluster = "cluster",
  weight = "weight",
  cfg = cfg,
  verbose = FALSE
)

res_complex
```

# Domain estimation

The `division` argument generates separate results for each observed category
of a division variable.

```{r domain-estimation}
res_domain <- svySE_calc(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  strata = "strata",
  cluster = "cluster",
  weight = "weight",
  division = "service",
  div_weight = NULL,
  cfg = cfg,
  verbose = FALSE
)

names(res_domain$results$ind_1$error)
```

The `TOTAL` element contains the complete estimate, while the remaining
elements contain the corresponding division results.

```{r inspect-domain}
res_domain$results$ind_1$error$S1
```

When a division requires a different expansion factor, it can be provided
through `div_weight`.

# Calculate simple indicator tables

`svySE_simple()` calculates unweighted frequencies and percentages without
requiring sampling weights, strata, or clusters.

```{r simple-tables}
res_simple <- svySE_simple(
  data = df,
  indicators = c("ind_1", "ind_2"),
  group_vars = "dept",
  group_labels = "Department",
  division = NULL,
  target = 1,
  valid_values = c(0, 1),
  pct_mult = 100,
  na_rm = TRUE,
  verbose = FALSE
)

res_simple
```

# Missing values in simple indicator tables

`svySE_simple()` provides explicit control over missing indicator values through
the `na_rm` argument.

When `na_rm = TRUE`, missing values are removed before frequencies and
percentages are calculated.

```{r simple-missing-values}
df_missing <- data.frame(
  dept = c("A", "A", "B", "B", "C", "C"),
  ind_1 = c(1, 0, NA, NA, 1, 1)
)

res_missing <- svySE_simple(
  data = df_missing,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  na_rm = TRUE,
  verbose = FALSE
)

res_missing$results$ind_1$simple$TOTAL
```

Department `B` is omitted because it does not contain any valid observations
for `ind_1`.

Missing observations are also excluded from the national denominator.

When `na_rm = FALSE`, missing indicator values are treated as a validation
condition and the calculation stops.

```{r simple-missing-error, error=TRUE}
svySE_simple(
  data = df_missing,
  indicators = "ind_1",
  group_vars = "dept",
  na_rm = FALSE,
  verbose = FALSE
)
```

This behavior makes the treatment of missing indicator values explicit and
prevents missing observations from being silently included in the denominator.

The returned object has class `"svySE_simple_result"`.

```{r simple-class}
class(res_simple)
```

A simple table can be inspected with:

```{r inspect-simple}
res_simple$results$ind_1$simple$TOTAL
```

The table contains:

| Column | Description |
|--------|-------------|
| `freq_0` | Frequency of non-target cases |
| `pct_0` | Percentage of non-target cases |
| `freq_1` | Frequency of target cases |
| `pct_1` | Percentage of target cases |
| `freq_total` | Total number of valid observations |
| `pct_total` | Total percentage |

The results generated by `svySE_simple()` describe the observed sample only.
Because no sampling weights or design variables are used, these percentages
should not be interpreted as population estimates.

# Simple tables by division

A division variable can also be used with `svySE_simple()`.

```{r simple-domain}
res_simple_domain <- svySE_simple(
  data = df,
  indicators = "ind_1",
  group_vars = "dept",
  group_labels = "Department",
  division = "service",
  target = 1,
  valid_values = c(0, 1),
  pct_mult = 100,
  verbose = FALSE
)

names(res_simple_domain$results$ind_1$simple)
```

# Export sampling errors to XLSX

The following example writes a sampling error workbook to a temporary file.

```{r export-error}
file_err <- tempfile(fileext = ".xlsx")

export_error <- svySE_xlsx(
  x = res_error,
  file_err = file_err,
  file_tab = NULL,
  cols_err = svySE_cols_err("full"),
  overwrite = TRUE
)

file.exists(file_err)
```

# Export simple tables to XLSX

```{r export-simple}
file_tab <- tempfile(fileext = ".xlsx")

export_simple <- svySE_xlsx(
  x = res_simple,
  file_err = NULL,
  file_tab = file_tab,
  cols_tab = NULL,
  overwrite = TRUE
)

file.exists(file_tab)
```

# Export multiple analyses

`svySE_xlsx()` can receive a named list containing several result objects.

```{r multiple-results}
results <- list(
  Main_errors = res_error,
  Domain_errors = res_domain,
  Main_simple = res_simple,
  Domain_simple = res_simple_domain
)
```

The function automatically identifies the class of each object:

- `"svySE_result"` objects are exported to `file_err`;
- `"svySE_simple_result"` objects are exported to `file_tab`.

```{r export-multiple}
multiple_err <- tempfile(fileext = ".xlsx")
multiple_tab <- tempfile(fileext = ".xlsx")

export_multiple <- svySE_xlsx(
  x = results,
  file_err = multiple_err,
  file_tab = multiple_tab,
  cols_err = svySE_cols_err("full"),
  cols_tab = NULL,
  overwrite = TRUE
)

file.exists(multiple_err)
file.exists(multiple_tab)
```

This is useful when several indicators are calculated from different datasets,
survey weights, domains, or analytical processes.

# Export selected results

The `select` argument restricts the export to chosen elements of a named list.

## Select sampling error results

```{r select-errors}
selected_err <- tempfile(fileext = ".xlsx")

svySE_xlsx(
  x = results,
  select = c("Main_errors", "Domain_errors"),
  file_err = selected_err,
  file_tab = NULL,
  overwrite = TRUE
)

file.exists(selected_err)
```

## Select simple results

```{r select-simple}
selected_tab <- tempfile(fileext = ".xlsx")

svySE_xlsx(
  x = results,
  select = c("Main_simple", "Domain_simple"),
  file_err = NULL,
  file_tab = selected_tab,
  overwrite = TRUE
)

file.exists(selected_tab)
```

# Customize exported columns

## Sampling error columns

The complete set of sampling error columns can be requested with:

```{r error-columns-full}
svySE_cols_err("full")
```

A custom selection can be defined with:

```{r error-columns-custom}
error_columns <- svySE_cols_err(
  type = "custom",
  cols = c(
    "est_pct",
    "se_pct",
    "ci_l_pct",
    "ci_u_pct",
    "cv",
    "deff",
    "n_unw"
  )
)

error_columns
```

Use the custom selection during export:

```{r export-custom-errors}
custom_err <- tempfile(fileext = ".xlsx")

svySE_xlsx(
  x = res_error,
  file_err = custom_err,
  file_tab = NULL,
  cols_err = error_columns,
  overwrite = TRUE
)

file.exists(custom_err)
```

## Simple table columns

```{r simple-columns-full}
svySE_cols_tab("full")
```

A custom selection can be defined with:

```{r simple-columns-custom}
simple_columns <- svySE_cols_tab(
  type = "custom",
  cols = c(
    "freq_1",
    "pct_1",
    "freq_total"
  )
)

simple_columns
```

Use the custom selection during export:

```{r export-custom-simple}
custom_tab <- tempfile(fileext = ".xlsx")

svySE_xlsx(
  x = res_simple,
  file_err = NULL,
  file_tab = custom_tab,
  cols_tab = simple_columns,
  overwrite = TRUE
)

file.exists(custom_tab)
```

# Main functions

| Function | Description |
|----------|-------------|
| `svySE_cfg()` | Configure sampling error estimation |
| `svySE_calc()` | Calculate weighted estimates and sampling errors |
| `svySE_simple()` | Calculate unweighted frequencies and percentages |
| `svySE_xlsx()` | Export one or multiple results to `.xlsx` files |
| `svySE_cols_err()` | Select sampling error columns |
| `svySE_cols_tab()` | Select simple table columns |

# Technical foundation

`svySE` uses established R packages for its internal workflow.

| Package | Role |
|---------|------|
| `survey` | Design-based estimation and variance calculation |
| `openxlsx` | Creation and formatting of `.xlsx` workbooks |
| `stats` | Statistical formulas, coefficients, and confidence intervals |
| `svySE` | High-level workflow for indicators, sampling errors, tables, and export |

`svySE` does not replace `survey`. It provides a structured interface for
repeated indicator production, quality measurement, and export.

# Summary

This vignette demonstrated how to:

1. configure sampling error estimation with `svySE_cfg()`;
2. calculate weighted estimates and sampling errors with `svySE_calc()`;
3. choose between Wald and logit (`xlogit`) confidence intervals for proportions;
4. work with weight-only, stratified, clustered, and combined designs;
5. calculate domain estimates;
6. generate unweighted descriptive tables with `svySE_simple()`;
7. export sampling errors and simple tables separately;
8. consolidate several analyses in one export process;
9. select specific results before exporting;
10. customize the columns included in `.xlsx` files.

The same workflow can be applied to official survey datasets by specifying the
corresponding indicators, grouping variables, weights, strata, and clusters.

## Exporting simple-table column profiles

When `cols_tab = NULL`, `svySE_xlsx()` exports all simple-table columns actually
available in each result object:

```{r, eval=FALSE}
svySE_xlsx(
  x = res_both,
  file_tab = "simple_all.xlsx",
  cols_tab = NULL
)
```

Predefined profiles are available:

```{r}
svySE_cols_tab("unweighted")
svySE_cols_tab("freq")
svySE_cols_tab("pct")
svySE_cols_tab("expanded")
svySE_cols_tab("expanded_freq")
svySE_cols_tab("expanded_pct")
svySE_cols_tab("counts")
svySE_cols_tab("percentages")
```

For example, export only sample and expanded counts:

```{r, eval=FALSE}
svySE_xlsx(
  x = res_both,
  file_tab = "simple_counts.xlsx",
  cols_tab = svySE_cols_tab("counts")
)
```

A fully custom selection is also supported:

```{r, eval=FALSE}
custom_columns <- svySE_cols_tab(
  type = "custom",
  cols = c(
    "freq_1",
    "pct_1",
    "exp_1",
    "exp_pct_1",
    "freq_total",
    "exp_total",
    "exp_pct_total"
  )
)

svySE_xlsx(
  x = res_both,
  file_tab = "simple_custom.xlsx",
  cols_tab = custom_columns
)
```

The requested columns must exist in the calculated object. For example,
expanded columns require `output = "weighted"` or `output = "both"`.