Estimate variance components

Estimate variance components for multistage sampling from frame or sample data, or infer a measure of homogeneity from a previous round's published design effect. Supports SRS and PPS selection of PSUs and, at three stages, of SSUs.

Usage

varcomp(x, ...)

# S3 method for class 'formula'
varcomp(
  x,
  ...,
  data = NULL,
  prob = NULL,
  strata = NULL,
  weights = NULL,
  estimator = c("expansion", "ratio")
)

# Default S3 method
varcomp(
  x,
  ...,
  stage_id = NULL,
  prob = NULL,
  strata = NULL,
  weights = NULL,
  estimator = c("expansion", "ratio"),
  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,
  estimator = c("expansion", "ratio")
)

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. For three stages, a list of two such entries gives the PSU and the SSU stage, each NULL for SRS (e.g., list(~psu_share, ~ssu_share)). The SSU entry has one value per observation, constant within each SSU, and sums to 1 across the SSUs of each PSU. It is supported on a complete frame only, not with weights.
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). For a three-stage sample, two columns (~w2 + w3, or a two-column matrix or data frame) give the stage-2 and stage-3 conditional weights separately, which the element-stage correction needs when many of each SSU's elements are sampled. NULL (default) computes frame components.
estimator: The estimator the components describe, "expansion" (default) for an estimated total over a known population size, or "ratio" for an estimated total over an estimated number of units, the usual estimator of a mean or a proportion. The two differ only when PSU sizes vary and PSUs are not drawn proportional to size. See Details. Not used on the back-out.
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. A published DEFT is squared first, and the published survey's weighting loss may need dividing out. 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: Take per PSU the published design effect refers to, counted in the units the estimate was computed on. For a person-level estimate it is the responding eligible persons per cluster, not the households per cluster. Give 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.
selection: How each sampled stage above the units is drawn, "srs" or "pps", named psu and, at three stages, ssu. NA for a backed-out icc, where the design effect does not say.
input: What the components were computed from: "frame", "weighted" (a sample with within-cluster weights), "deff" or "se" (a published design effect or standard error).
estimator: The estimator the components describe, "expansion" or "ratio". NA for a backed-out icc.

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, and 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.

The element-stage noise in a 3-stage sample depends on the share of each SSU's elements that was sampled, $f_3 = n_{ij} / N_{ij}$. Only the stage-3 weight sums to $N_{ij}$, so a single within-cluster weight column reads the fraction as $f_2 f_3$ and removes too much noise from the SSU totals, which understates varw_psu and icc_ssu. The difference is small when $f_3$ is small, the usual case of a few households listed from many. When it is not, pass the stage-2 and stage-3 weights as two columns. At $f_2 = 0.2$ and $f_3 = 0.5$ one column understates varw_psu by about a fifth.

A sample passed without weights is read as the frame, and that overstates icc. Each sampled PSU's mean carries the noise of its own subsample, about $(1 - f_2) S_w^2 / m$ for a take of $m$ from a PSU with within-PSU variance $S_w^2$, and the between-PSU component absorbs it, so planning with it counts the within-PSU variance twice. The weights correction is what removes that term. Neither path estimates the sampling variance of the sample's own estimate, whose estimator combines the two components differently.

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. At three stages, ids = ~psu + ssu, a third probability column is the element stage and supplies the stage-3 weight on its own. Without one, every element of a sampled SSU is taken, so the SSU totals carry no element-stage noise. 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, contributing no between-PSU variance and needing to be treated as separate strata, with prob renormalized over the remaining PSUs.

Planning from a complete PSU register, as n_alloc() does with psu, is the exception. The certainty PSUs are not known until the allocation is, so use every PSU of the stratum, with prob equal to each PSU's share of the stratum's units. The icc returned is then the size-weighted homogeneity that the allocation's remainder model assumes. The SRS value measures the spread of PSU totals, which carries the spread of PSU sizes, so with unequal sizes it is not the homogeneity that model needs.

For two-stage designs, icc represents the measure of homogeneity $\delta = V_b / (V_b + V_w)$, written $\delta$ in Valliant, Dever, and Kreuter (2018, Ch. 9), in the multistage sampling variance decomposition. It lies in $[0, 1]$ when the components are nonnegative and their sum is positive. It is related to conventional intraclass correlation measures, including Kish's rate of homogeneity (roh), but is not generally numerically identical to them. Survey ICC estimates such as roh can be negative. A random-intercept ICC formed from nonnegative variance components cannot. When using an ICC from another source, check its definition and the associated design-effect formula. Here the two-stage approximation is $k[1 + (b - 1)\delta]$, with var_ratio supplying $k$. The three-stage component definitions are given below.

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_, not being identified. Functions that need the unit relvariance, such as n_cluster() with a cv target, say so and ask for it directly. Inverting a published design effect can return an icc outside $[0, 1]$. The value is retained with a warning, but the cluster-planning functions do not accept it.

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$, so 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 overstates a positive $\delta$ by roughly the same relative amount that $b^*$ exceeds the nominal figure, about 9 percent at a 30 percent coefficient of variation (0.0421 against 0.0385 at a design effect of 1.8 and a take of 20). 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.

A published DEFT and its weighting loss

deff is a ratio of variances. Some survey programs publish DEFT, the ratio of standard errors, which is its square root. The DHS sampling-error appendix lists it under the heading "Design effect (DEFT)", so square the published value before passing it. MICS tables give both columns.

A published design effect also carries the weighting loss of the survey that produced it, which the identity above attributes to $\delta$. With that survey's final weights w, Kish's factor from design_effect() divides it out:

kish <- design_effect(weights = w)
vc <- varcomp(deff = deft^2 / kish, n_per_psu = b)

This is the decomposition Kalton, Brick and Lê (2005) use to carry a design effect from one survey to the planning of another, with the new design's weighting factor multiplied back in. Gabler, Haeder and Lahiri (1999) justify the product of the two factors under a model, as a conservative value for the design effect. The division assumes the weights are unrelated to the outcome, as Kish's decomposition does. Passing the factor as var_ratio returns the same icc, but the factor then stays in the result and multiplies every design effect formed from it.

Whether to divide depends on how the new plan is built. n_alloc() over the strata forms the variance of a disproportionate allocation from the allocation itself, so the old factor is divided out, or it would count twice. A scalar deff in n_prop() or n_mean() carries no allocation. For a new design as disproportionate as the old one, multiply its expected weighting factor back in, deff = design_effect(vc, n_per_psu = b_new) * kish_new, or the sample comes out too small. Weight variation within strata, from nonresponse adjustment for instance, recurs in the new survey and n_alloc() does not model it. A deff column in its frame can carry it.

A published standard error carries the same loss. To remove it, form the design effect and divide, deff = se^2 / (p * (1 - p) / n) / kish, rather than passing se.

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.

These are the components of the expansion estimator, which divides an estimated total by a known population size. A mean or a proportion is usually estimated as a ratio instead, the estimated total over the estimated number of units, and with PSUs of unequal size drawn by SRS its variance can be much smaller. The PSU totals of the expansion estimator vary with the PSU sizes as well as with their content, and the ratio estimator removes the part due to size. estimator = "ratio" returns its components, the same expressions with each total $t_i$ replaced by $t_i - N_i\bar y$ for the overall mean $\bar y$, at every stage. The within-PSU variances and unit_relvar are unchanged. With PSUs of equal size, or a first stage drawn with probability proportional to size, the two estimators give the same components. Use "ratio" when planning a mean or a proportion from PSUs of unequal size drawn with equal probability, and the default when planning a total.

For three stages, varb uses the between-PSU expression above that matches the first stage, the SRS one without prob and the PPS one with it. varw_psu scales the variance of SSU totals within each PSU by the inverse draw probability $1/p_i$, which is $M$ without prob, 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, and var_ratio_ssu depends on how SSUs are drawn. By default SSUs are taken by simple random sampling within each PSU and their totals are expanded by the number of SSUs. When SSU sizes vary, those totals vary with size as well as with content, so var_ratio_ssu rises above the var_ratio_psu * (1 - icc_psu) that planning functions derive when it is not supplied, and the more so the smaller the unit relvariance. On a frame of census tracts and block groups this ranged from 1.2 to 3.2 against derived values near 1. A second-stage prob (one-draw probabilities of the SSUs within their PSU, usually their size share) describes selection proportional to size instead. The size effect then cancels, and var_ratio_ssu comes close to the derived value. Use the prob that matches the design being planned. Passing a svyplan_varcomp to a planning function uses the estimated pair as given.

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.

icc and the intraclass correlation

icc is the between-PSU share of the two components, so it is never negative. The intraclass correlation of single-stage cluster sampling, the average correlation between units of the same cluster, can be negative, and a cluster sample is then more efficient than a simple random sample of the same size. In this parameterization that case shows as a small positive icc, not a negative one. Under the stage-by-stage correction (fpc = "stage" in n_cluster()), a design whose icc <= 1 / (M + 1), for M units per PSU, is most efficient taking whole PSUs. For PSUs of equal size drawn by simple random sampling, and frame components, whole PSUs beat a simple random sample of the same size exactly when icc < 1 / (M + 1). A negative correlation published elsewhere therefore cannot be passed as icc. Estimate the components from the frame instead.

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.

Kalton, G., Brick, J. M., and Lê, T. (2005). Estimating components of design effects for use in sample design. In Household Sample Surveys in Developing and Transition Countries, Studies in Methods, Series F No. 96, Chapter VI. United Nations.

Gabler, S., Haeder, S., and Lahiri, P. (1999). A model based justification of Kish's formula for design effects for weighting and clustering. Survey Methodology, 25(1), 105-106.

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, SRS PSUs, frame)
#> 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, PPS PSUs, weighted sample)
#> 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, SRS PSUs, frame)
#> 
#>  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, SRS PSUs, SRS SSUs, frame)
#> varb = 0.0016, varw_psu = 0.0111, varw_ssu = 0.0442
#> icc = 0.0347, 0.2002
#> var_ratio = 1.0500, 1.2670
#> Unit relvariance = 0.0436

# A 3-stage sample: 2 of 5 villages per district, 2 of 4 units per
# village. Half of each village is sampled, so the stage-2 and stage-3
# weights go in as two columns
villages <- unlist(lapply(split(unique(frame3$village),
                                rep(1:20, each = 5))[1:6],
                          sample, size = 2))
units <- unlist(lapply(split(seq_len(400), frame3$village)[villages],
                       sample, size = 2))
samp3 <- frame3[units, ]
samp3$w2 <- 5 / 2
samp3$w3 <- 4 / 2
varcomp(income ~ district/village, data = samp3, weights = ~ w2 + w3)
#> Variance components (3-stage, SRS PSUs, SRS SSUs, weighted sample)
#> varb = 0.0193, varw_psu = 0.0412, varw_ssu = 0.0276
#> icc = 0.2769, 0.5990
#> var_ratio = 1.1770, 1.1602
#> Unit relvariance = 0.0593

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

# Villages of unequal size. Drawn by SRS, their totals vary with size
# and var_ratio_ssu is large. Drawn proportional to size (each
# village's share of its district), it is close to the derived value.
set.seed(42)
sizes <- rpois(100, 8) + 2
frame4 <- data.frame(
  district = rep(rep(1:20, each = 5), sizes),
  village = rep(1:100, sizes),
  income = rnorm(sum(sizes), 50000, 10000)
)
n_district <- ave(frame4$income, frame4$district, FUN = length)
frame4$district_share <- n_district / nrow(frame4)
frame4$village_share <- ave(frame4$income, frame4$village, FUN = length) /
  n_district
varcomp(income ~ district/village, data = frame4,
        prob = ~district_share)$var_ratio
#> var_ratio_psu var_ratio_ssu 
#>      1.018853      3.540807 
vc4 <- varcomp(income ~ district/village, data = frame4,
               prob = list(~district_share, ~village_share))
vc4$var_ratio
#> var_ratio_psu var_ratio_ssu 
#>      1.018853      1.080727 
vc4$var_ratio[[1]] * (1 - vc4$icc[[1]])
#> [1] 1.002051

# Vector (list) interface
varcomp(frame3$income,
        stage_id = list(frame3$district, frame3$village))
#> Variance components (3-stage, SRS PSUs, SRS SSUs, frame)
#> varb = 0.0016, varw_psu = 0.0111, varw_ssu = 0.0442
#> icc = 0.0347, 0.2002
#> var_ratio = 1.0500, 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

# A DHS appendix publishes DEFT, the ratio of standard errors
deft <- 1.582
varcomp(deff = deft^2, n_per_psu = 12.3)
#> Variance components (2-stage, from a design effect)
#> icc = 0.1330
#> var_ratio = 1.0000
#> deff = 2.5027 at n_per_psu = 12.3
#> varb, varw and unit_relvar are not identified by a design effect

# Divide out the weighting loss of the survey that published it
w_old <- rep(c(1, 4), c(300, 100))
kish <- design_effect(weights = w_old)
vc_dhs <- varcomp(deff = deft^2 / kish, n_per_psu = 12.3)
vc_dhs
#> Variance components (2-stage, from a design effect)
#> icc = 0.0543
#> var_ratio = 1.0000
#> deff = 1.6136 at n_per_psu = 12.3
#> varb, varw and unit_relvar are not identified by a design effect

# A scalar deff for an equally disproportionate design puts it back
n_prop(p = 0.28, cv = 0.05,
       deff = design_effect(vc_dhs, n_per_psu = 15) * kish)
#> Sample size for proportion (wald)
#> n = 2809 (p = 0.28, cv = 0.050, deff = 2.73)
#> expected cases = 786.3