Skip to contents

Design-weighted estimates from a probability sample

Source: R/core-stats-estimate.R
hr_ht_estimate.Rd

Estimates a population total, mean or ratio from units sampled with known inclusion probabilities pi, such as a reference that judged every kept record and a 10% sample of the excluded ones.

Usage

hr_ht_estimate(
  y,
  pi,
  strata = NULL,
  estimand = c("mean", "total", "ratio"),
  den = NULL,
  yhat = NULL,
  yhat_pop = NULL,
  B = 2000,
  conf = 0.95,
  seed = 1
)

Arguments

y

Values for the sampled units (numeric or logical), no NA.

pi

Inclusion probability of each sampled unit, in (0, 1].

strata

Optional stratum of each sampled unit; the bootstrap resamples within strata.

estimand

"mean", "total" or "ratio".

den

For "ratio", the denominator variable on the sampled units.

yhat, yhat_pop

Predictions on the sampled units and on every unit in the population, for the difference estimator. Give both or neither. Not used with "ratio".

B

Number of bootstrap resamples; 0 skips the interval.

conf

Confidence level.

seed

Seed for the bootstrap.

Value

A one-row data frame: estimand, method ("horvitz_thompson", "hajek", "ratio" or "difference"), estimate, se, lower, upper, n (sampled units), N_hat (sum(1 / pi), the estimated population size), conf, B and note.

Details

  • "total": the Horvitz-Thompson estimator, sum(y / pi).

  • "mean": the Hajek estimator, sum(y / pi) / sum(1 / pi), which divides by the estimated population size and is steadier than the HT total over a known N when the sample size varies.

  • "ratio": the ratio of two HT totals, sum(y / pi) / sum(den / pi). With y = "eligible and kept" and den = "eligible", this is design-weighted recall.

With yhat (predictions for the sampled units) and yhat_pop (predictions for every unit in the population) the difference estimator is used instead: total = sum(yhat_pop) + sum((y - yhat) / pi), and the mean is that total over length(yhat_pop). It stays unbiased whatever the quality of the predictions, and its variance shrinks as they improve.

The interval is a percentile bootstrap with the finite population correction: the rescaling bootstrap of Rao and Wu (1988) with n_h - 1 draws per stratum. In each of B resamples (under seed, with the caller's random state put back) each stratum h draws n_h - 1 of its n_h sampled units with replacement, and unit i gets the weight multiplier 1 - lambda_h + lambda_h * r_i * n_h / (n_h - 1), where r_i is the number of times it was drawn and lambda_h = sqrt(1 - f_h). The sampling fraction f_h is n_h / sum(1 / pi) over the stratum, which is pi itself when pi is constant within the stratum. For an HT total under stratified simple random sampling the replicate variance then has expectation sum_h N_h^2 (1 - f_h) s_h^2 / n_h, the design variance. A stratum judged in full (pi = 1, so f_h = 1) is a census and adds no variance whatever its values. se is the standard deviation of the replicates. A stratum with one sampled unit also adds no variance, and the result notes it. Give strata that match the design, with one pi per stratum, for the correction to be exact.

References

Horvitz, D. G. and Thompson, D. J. (1952). A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association 47(260), 663-685.

Hajek, J. (1971). Comment on a paper by D. Basu. In Godambe, V. P. and Sprott, D. A. (eds), Foundations of Statistical Inference, 236.

Rao, J. N. K. and Wu, C. F. J. (1988). Resampling inference with complex survey data. Journal of the American Statistical Association 83(401), 231-241.

Sarndal, C.-E., Swensson, B. and Wretman, J. (1992). Model Assisted Survey Sampling. Springer. (The difference estimator, section 6.3.)

See also

Other statistics: hr_agreement(), hr_alt_test(), hr_ppi(), hr_prop_ci()

Examples

# Every kept record judged (pi = 1), one in ten excluded records judged.
eligible <- c(rep(1, 18), rep(0, 2), rep(1, 1), rep(0, 29))
kept <- c(rep(1, 20), rep(0, 30))
pi <- ifelse(kept == 1, 1, 0.1)

# Design-weighted recall: eligible and kept over all eligible.
hr_ht_estimate(eligible * kept, pi, strata = kept, estimand = "ratio",
               den = eligible, B = 500)
#>   estimand method  estimate        se     lower    upper  n N_hat conf   B note
#> 1    ratio  ratio 0.6428571 0.2200034 0.3753514 0.972281 50   320 0.95 500 <NA>

# Estimated number of eligible records in the population.
hr_ht_estimate(eligible, pi, strata = kept, estimand = "total", B = 500)
#>   estimand           method estimate       se    lower    upper  n N_hat conf
#> 1    total horvitz_thompson       28 9.611748 18.51317 47.95506 50   320 0.95
#>     B note
#> 1 500 <NA>