Estimate variance components

Estimate between- and within-stage variance components using nested ANOVA decomposition, from frame data, from a sample, or from a previous round's published design effect. 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,
  deff = NULL,
  se = NULL,
  p = NULL,
  var = NULL,
  n = NULL,
  n_per_psu = NULL,
  n_per_ssu = NULL,
  cv_take = NULL,
  var_ratio = NULL
)

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

Arguments

x: A formula, numeric vector, or survey design object (see Details). Leave unset to back an icc out of a published design effect.
...: 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).
deff: A published design effect to back an icc out of, with x unset. Requires n_per_psu. See Details.
se: A published standard error, an alternative to deff for the same back-out. Requires n and one of p or var.
p: Proportion the published se refers to.
var: Unit variance the published se refers to, for a mean.
n: Sample size the published se was computed on.
n_per_psu: Units per PSU the published design effect refers to: the size-weighted average take, a vector of realized takes, or a nominal take with cv_take. See Details.
n_per_ssu: Units per SSU. Not accepted on the back-out, which a scalar design effect does not identify.
cv_take: Coefficient of variation of the realized takes, when n_per_psu is their nominal mean.
var_ratio: Ratio of the clustered component's unit variance to the analysis variable's, 1 when unset. See design_effect().

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.
icc: Measure of homogeneity. Length 1 for 2-stage, length 2 (icc_psu, icc_ssu) for 3-stage.
var_ratio: Ratio parameter(s), same length as icc (var_ratio_psu, var_ratio_ssu for 3-stage).
unit_relvar: Unit relvariance (scalar).
stages: Number of stages (2 or 3).
strata: Per-stratum component table when strata is supplied, otherwise NULL (see Details).
source: "data" when the components were estimated, "deff" when icc was backed out of a design effect.
params: The take that identifies a backed-out icc, n_per_psu alongside the nominal n_per_psu_nominal. Empty on the estimated path, where the take is a later design choice.

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. A backed-out result adds n_per_psu and source to that row, and carries NA in the three fields a design effect does not identify.

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, first-stage strata, and within-cluster weights are extracted from the design object; strata and weights override what is extracted. Requires the survey package. For a PPS first stage, also pass the PSU selection probabilities via prob. Which designs can supply their own weights is set out under Details.

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 icc_psu) is conservative: biased upward, never downward.

What a survey.design can supply on its own

The survey.design method needs those same within-cluster weights, and whether it can recover them depends on how the design was built rather than on how deep its ids goes. A design built with per-stage probabilities, svydesign(ids = ~psu, probs = ~p1 + p2), keeps the stages apart, and the columns after the first multiply to exactly the weight required. A design built with weights= stores one combined weight per element however many stages ids names, and the stage-1 factor cannot be divided back out of it.

Unit weights describe a frame and pass straight through. Any other single-column weight is refused rather than guessed at, because a design declared as ids = ~psu carrying a full design weight is equally consistent with a census inside each selected cluster and with a two-stage design that folded its second stage into the weight, and those give materially different components. Supply the conditional weights through weights= in that case, as a one-sided formula naming a design variable or as a numeric vector.

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 icc is the design-based measure of homogeneity $\delta = V_b / (V_b + V_w)$, written $\delta$ in Valliant, Dever, and Kreuter (2018, Ch. 9). Unlike the traditional ANOVA intraclass correlation coefficient, icc is constrained to $[0, 1]$ and should not be compared directly to mixed-model ICCs (e.g. from lme4) which can be negative.

From a published design effect

The most widely available source of an icc is a previous round's report, which publishes a design effect and an average cluster take rather than the microdata the other interfaces need. Leaving x unset inverts the two-stage clustering identity used by design_effect(), $$D = k(1 + \delta(b - 1)) \quad\Longrightarrow\quad \delta = (D/k - 1)/(b - 1),$$ and returns the icc in a svyplan_varcomp the cluster functions accept:

vc <- varcomp(deff = 1.8, n_per_psu = 20)
n_cluster(stage_cost = c(500, 50), icc = vc, budget = 1e5)

A scalar design effect fixes one product and nothing more, so varb, varw and unit_relvar come back NA_real_: they are not identified. Functions that need the unit relvariance, such as n_cluster() with a cv target, say so and ask for it directly.

Where the report gives a standard error instead of a design effect, se with the quantity it refers to forms the design effect against the simple random sample of the same size, $D = s^2 / (pq/n)$ for a proportion and $D = s^2 / (S^2/n)$ for a mean.

Which take identifies the icc

The $b$ in the identity is the size-weighted average number of sampled units per PSU, $$b^* = \sum_i b_i^2 / \sum_i b_i = \bar b (1 + cv_b^2),$$ with $cv_b$ the population coefficient of variation of the realized takes. A planned take is a constant, so $cv_b = 0$ and $b^* = b$: the forward direction in design_effect() never meets the distinction. A published design effect was achieved with takes that varied, and since $\delta$ moves inversely with $b - 1$, using the nominal take understates it by roughly the same relative amount that $b^*$ exceeds the nominal figure, about 9 percent at a 30 percent coefficient of variation. Report tables carry the nominal figure, which is the one a reader reaches for, so supply the realized takes or their coefficient of variation where they are known:

varcomp(deff = 1.8, n_per_psu = 20)                  # constant take
varcomp(deff = 1.8, n_per_psu = takes)               # realized takes
varcomp(deff = 1.8, n_per_psu = 20, cv_take = 0.3)   # summary form

The take is stored in $params as provenance, since the same icc means different things at a take of 12 and of 30, and design_effect() does not read it back: the reason to back an icc out at all is to re-plan at a different take, and defaulting to the old one would return the old design's number.

Component estimands

For two-stage SRS, let $M$ be the number of PSUs, $N_i$ the ultimate-unit count in PSU $i$, $t_i$ its outcome total, and $S_i^2$ its within-PSU variance. With $t_U=M\bar t$, the complete-frame components returned are $$V_b=s^2(t_i)/\bar t^2, \qquad V_w=M\sum_iN_i^2S_i^2/t_U^2.$$ For PPS with one-draw probabilities $p_i$, they are $$V_b=\sum_i p_i(t_i/p_i-t_U)^2/t_U^2, \qquad V_w=\sum_iN_i^2S_i^2/(p_it_U^2).$$ In both cases unit_relvar is the ultimate-unit variance divided by the squared ultimate-unit mean, icc = varb / (varb + varw), and var_ratio = (varb + varw) / unit_relvar.

For three stages, varb uses the PPS between-PSU expression above. varw_psu applies the same PPS scaling to the variance of SSU totals within each PSU, and varw_ssu aggregates the within-SSU ultimate-unit variances. The returned PSU homogeneity compares the between-PSU component with the element-level variance aggregated within PSUs; the SSU homogeneity is varw_psu / (varw_psu + varw_ssu). This distinction is why the first icc is not generally varb / sum(c(varb, varw)) for a three-stage result.

The two var_ratio values are estimated from their own stage decompositions rather than by imposing the identity var_ratio_ssu = var_ratio_psu * (1 - icc_psu) that the planning formula uses (see design_effect()). On small clusters the two can differ by several percent. Passing a svyplan_varcomp to a planning function uses the estimated pair as given; omitting var_ratio_ssu there applies the identity instead.

When weights are supplied, $N_i$, totals, means, and variances in these expressions are their weighted estimates. The code subtracts the stated SRS-without-replacement estimation variance of estimated cluster totals and truncates negative corrected components to zero. These corrections are planning approximations, not exact variance estimators for arbitrary informative multistage samples.

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, icc, ...) are not filled. The table's columns (sd, mean, icc_psu, var_ratio_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, and for the back-out from a published design effect (leave x unset).
varcomp(survey.design): Method for survey design objects. Pass a one-sided formula (e.g., ~y) to specify the outcome variable. Cluster structure, first-stage strata, and within-cluster weights are extracted from the design; strata and weights override what is extracted.

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
#> icc = 0.1252
#> var_ratio = 1.0979
#> Unit relvariance = 0.0365
as.data.frame(vc2)
#>   stages        varb       varw       icc var_ratio unit_relvar
#> 1      2 0.005016206 0.03504851 0.1252026  1.097919  0.03649151

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

# 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
#> icc = 0.1115
#> var_ratio = 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       icc_psu var_ratio_psu varb   varw   unit_relvar
#>  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
#> icc = 0.0330, 0.2002
#> var_ratio = 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
#> icc = 0.0330, 0.2002
#> var_ratio = 1.0482, 1.2670
#> Unit relvariance = 0.0436

# Back an icc out of a previous round's published design effect
vc <- varcomp(deff = 1.8, n_per_psu = 20)
vc
#> Variance components (2-stage, from a design effect)
#> icc = 0.0421
#> var_ratio = 1.0000
#> deff = 1.8000 at n_per_psu = 20
#> varb, varw and unit_relvar are not identified by a design effect

# Re-plan at a different take: the reason to back the icc out
design_effect(vc, n_per_psu = 12)
#> Planning design effect: 1.4632

# Realized takes varied around a nominal 20
varcomp(deff = 1.8, n_per_psu = 20, cv_take = 0.3)
#> Variance components (2-stage, from a design effect)
#> icc = 0.0385
#> var_ratio = 1.0000
#> deff = 1.8000 at n_per_psu = 21.8 (nominal 20, size-weighted)
#> varb, varw and unit_relvar are not identified by a design effect

# From a published standard error instead
varcomp(se = 0.021, p = 0.30, n = 1200, n_per_psu = 20)
#> Variance components (2-stage, from a design effect)
#> icc = 0.0800
#> var_ratio = 1.0000
#> deff = 2.5200 at n_per_psu = 20
#> varb, varw and unit_relvar are not identified by a design effect