Estimate Variance Components

Estimate between- and within-stage variance components using nested ANOVA decomposition. Supports SRS and PPS first-stage designs.

Usage

varcomp(x, ...)

# S3 method for class 'formula'
varcomp(x, ..., data = NULL, prob = NULL, strata = NULL, weights = NULL)

# Default S3 method
varcomp(x, ..., stage_id = NULL, prob = NULL, strata = NULL, weights = NULL)

# S3 method for class 'survey.design'
varcomp(x, ..., prob = NULL, strata = NULL)

Arguments

x: A formula, numeric vector, or survey design object (see Details).
...: Additional arguments passed to methods. Unused arguments are rejected.
data: A data frame (required for formula interface).
prob: First-stage selection probabilities (PPS). A one-sided formula (e.g., ~pp) when using the formula interface, or a numeric vector: either one value per observation (constant within each PSU) or one value per PSU. A one-per-PSU vector is matched by name when named. Unnamed values are taken in sorted order of the unique PSU identifiers. Values must be strictly between 0 and 1 and sum to 1 across PSUs (tolerance 1e-3). NULL (default) assumes SRS.
strata: Optional stratification: a one-sided formula (formula and survey.design interfaces) or a vector (default interface). Components are then estimated per stratum. See Details.
weights: Optional sampling weights for estimating components from a sample rather than a frame: a one-sided formula (e.g., ~w) when using the formula interface, or a numeric vector with one positive value per observation. Within each cluster the summed weights must estimate the cluster population size. For a multi-stage sample use the within-cluster weights, not the full design weight (see Details). NULL (default) computes frame components.
stage_id: A list of stage-ID vectors (required for vector interface). Length determines the number of stage boundaries (stages - 1).

Value

A svyplan_varcomp object with components:

varb: Between-PSU variance (scalar).
varw: Within-PSU variance. Scalar for 2-stage, length-2 vector (varw_psu, varw_ssu) for 3-stage.
delta: Measure of homogeneity. Length 1 for 2-stage, length 2 (delta_psu, delta_ssu) for 3-stage.
k: Ratio parameter(s), same length as delta (k_psu, k_ssu for 3-stage).
rel_var: Unit relvariance (scalar).
stages: Number of stages (2 or 3).
strata: Per-stratum component table when strata is supplied, otherwise NULL (see Details).

The result can be exported with as.data.frame(). Unstratified results become a one-row two- or three-stage component table. Stratified results return the table stored in $strata.

Details

The interface is determined by the class of x:

Formula: varcomp(income ~ district/village, data = frame). The LHS is the analysis variable, RHS terms express nesting using / or %in% (see formula): y ~ psu/ssu or equivalently y ~ ssu %in% psu (SSUs nested within PSUs). With / the outermost stage comes first, whereas with %in% the innermost comes first. Be careful not to reverse the order.
Numeric vector: varcomp(y, stage_id = list(cluster_ids)).
survey.design: varcomp(design, ~y). Cluster structure and design weights are extracted from the design object. Requires the survey package. For a PPS first stage, also pass the PSU selection probabilities via prob.

The formula and vector interfaces compute frame components: data is treated as a complete population. To estimate components from a sample instead, supply weights (formula and vector interfaces) or use the survey.design method. Weights are treated as inverse inclusion probabilities: cluster sizes and totals are estimated by summed weights, and the estimation variance of the weighted cluster totals is subtracted from the between-stage variance terms, so unequal-probability samples from a previous round give approximately design-unbiased components. Unit weights recover the frame formulas exactly. The weight scale matters: within each cluster the summed weights should estimate the cluster population size, since the implied sampling fraction drives the correction. For a multi-stage sample this means the within-cluster weights (the product of the stage-2 and later weights), not the full design weight, whose stage-1 factor would overstate every cluster size. The correction assumes noninformative (SRS-like) subsampling within clusters. Informative within-cluster sampling remains approximate. In the 3-stage weighted case the between-PSU correction removes element-stage estimation noise but not SSU-stage subsampling noise, so when few SSUs are sampled per PSU the between-PSU component (and delta_psu) is conservative: biased upward, never downward.

When prob is NULL, SRS first-stage is assumed. When provided, PPS variance estimation is used. prob must sum to 1 across the PSUs present in the data. On a complete frame these are the one-draw selection probabilities themselves. In a sample containing only some of the PSUs, renormalize over the sampled PSUs: prob = pp[sampled] / sum(pp[sampled]). For a fixed-size PPS design the stage-1 inclusion probabilities are proportional to the one-draw probabilities, so shares derived from inverse stage-1 weights, pi1 / sum(pi1), are identical and can be used when the frame is not at hand. Certainty (take-all) PSUs have no place in the PPS path: they contribute no between-PSU variance and should be treated as separate strata, with prob renormalized over the remaining PSUs.

The returned delta is the measure of homogeneity $\delta = V_b / (V_b + V_w)$ following Valliant, Dever, and Kreuter (2018, Ch. 9). Unlike the traditional ANOVA intraclass correlation coefficient, delta is constrained to $[0, 1]$ and should not be compared directly to mixed-model ICCs (e.g. from lme4) which can be negative.

Clusters containing a single observation have undefined within-cluster variance. In this case, the within-cluster variance is imputed as the mean variance of the remaining clusters.

With strata, components are estimated separately within each stratum and returned as a per-stratum table in $strata (also via as.data.frame()). The pooled fields (varb, delta, ...) are not filled. The table's columns (sd, mean, delta_psu, k_psu) match the n_alloc() frame contract, so after adding stratum N it feeds a stratified two-stage allocation directly. When prob is combined with strata, supply one value per observation, summing to 1 within each stratum.

Methods (by class)

varcomp(formula): Method for formula interface.
varcomp(default): Default method for numeric vectors.
varcomp(survey.design): Method for survey design objects. Pass a one-sided formula (e.g., ~y) to specify the outcome variable. Cluster structure and design weights are extracted from the design.

References

Valliant, R., Dever, J. A., and Kreuter, F. (2018). Practical Tools for Designing and Weighting Survey Samples (2nd ed.). Springer. Ch. 9.

Hansen, M. H., Hurwitz, W. N., and Madow, W. G. (1953). Sample Survey Methods and Theory (Vol. I). Wiley.

Examples

# 2-stage SRS using formula (PSU = district)
set.seed(314)
frame2 <- data.frame(
  income = rnorm(200, 50000, 10000),
  district = rep(1:20, each = 10)
)
vc2 <- varcomp(income ~ district, data = frame2)
vc2
#> Variance components (2-stage)
#> varb = 0.0050, varw = 0.0350
#> delta = 0.1252
#> k = 1.0979
#> Unit relvariance = 0.0365
as.data.frame(vc2)
#>   stages        varb       varw     delta        k    rel_var
#> 1      2 0.005016206 0.03504851 0.1252026 1.097919 0.03649151

# Feed into n_cluster
n_cluster(stage_cost = c(500, 50), delta = vc2, budget = 100000)
#> Optimal 2-stage allocation
#> field design: n_psu = 111 | psu_size = 8 -> total n = 888
#> cv = 0.0092, cost = 99900
#> continuous optimum: n_psu = 108.9392 | psu_size = 8.35886 (cv = 0.0092, cost = 100000)

# Estimate components from a two-stage sample: within-cluster
# weights (summing to each cluster's population size) and, for a
# PPS first stage, one-draw probabilities renormalized over the
# sampled clusters
sampled <- unlist(lapply(split(seq_len(200), frame2$district)[1:8],
                         sample, size = 4))
samp <- frame2[sampled, ]
samp$w <- 10 / 4
pp <- rep(1 / 20, 8)
varcomp(income ~ district, data = samp, weights = ~w,
        prob = pp / sum(pp))
#> Variance components (2-stage)
#> varb = 0.0054, varw = 0.0428
#> delta = 0.1115
#> k = 1.0830
#> Unit relvariance = 0.0445

# Per-stratum components for a stratified two-stage plan
set.seed(2718)
frame_s <- data.frame(
  income = rnorm(400, 50000, 10000),
  district = rep(1:40, each = 10),
  region = rep(c("North", "South"), each = 200)
)
varcomp(income ~ district, data = frame_s, strata = ~region)
#> Variance components (2-stage, 2 strata)
#>  stratum sd         mean       delta_psu k_psu  varb   varw   rel_var
#>  North   10154.3236 49512.7588 0.0967    1.0997 0.0045 0.0418 0.0421 
#>  South   9553.0569  49952.3385 0.0779    1.1008 0.0031 0.0371 0.0366 

# 3-stage SRS using formula: villages nested within districts
# "/" expresses nesting (outermost stage first, see ?formula)
set.seed(1618)
frame3 <- data.frame(
  income = rnorm(400, 50000, 10000),
  district = rep(1:20, each = 20),
  village = rep(1:100, each = 4)
)
vc3 <- varcomp(income ~ district/village, data = frame3)
vc3
#> Variance components (3-stage)
#> varb = 0.0015, varw_psu = 0.0111, varw_ssu = 0.0442
#> delta = 0.0330, 0.2002
#> k = 1.0482, 1.2670
#> Unit relvariance = 0.0436

# 3-stage PPS (explicit first-stage probabilities)
frame3$pp <- rep(1 / 20, 400)
vc3_pps <- varcomp(income ~ district/village, data = frame3, prob = ~pp)

# Vector (list) interface
varcomp(frame3$income,
        stage_id = list(frame3$district, frame3$village))
#> Variance components (3-stage)
#> varb = 0.0015, varw_psu = 0.0111, varw_ssu = 0.0442
#> delta = 0.0330, 0.2002
#> k = 1.0482, 1.2670
#> Unit relvariance = 0.0436