HTvar: Variance of the Horvitz-Thompson estimator

In rhobis/jipApprox: Approximate Inclusion Probabilities for Survey Sampling

Variance of the Horvitz-Thompson estimator

Description

Compute or estimate the variance of the Horvitz-Thompson total estimator by the Horvitz-Thompson or Sen-Yates-Grundy variance estimators.

Usage

HTvar(y, pikl, sample = TRUE, method = "HT")

Arguments

y	numeric vector representing the variable of interest
pikl	matrix of second-order (joint) inclusion probabilities; the diagonal must contain the first-order inclusion probabilities.
sample	logical value indicating whether sample or population values are provided. If sample=TRUE, the function returns a sample estimate of the variance, while if sample=FALSE, the Variance is computed over all population units. Default is TRUE.
method	string, indicating if the Horvitz-Thompson ("HT") or the Sen-Yates-Grundy ("SYG") estimator should be computed.

Details

The Horvitz-Thompson variance is defined as

\sum_{i\in U}\sum_{j \in U} \frac{(\pi_{ij} - \pi_i\pi_j)}{\pi_i\pi_j} y_i y_j

which is estimated by

\sum_{i\in U}\sum_{j \in U} \frac{(\pi_{ij} - \pi_i\pi_j)}{\pi_i\pi_j\pi_{ij}} y_i y_j

The Sen-Yates-Grundy variance is obtained from the Horvitz-Thompson variance by conditioning on the sample size n, and is therefore only appliable to fixed size sampling designs:

\sum_{i\in U}\sum_{j > i} (\pi_i\pi_j - \pi_{ij}) \Biggl(\frac{y_i}{\pi_i} - \frac{y_j}{\pi_j} \Biggr)^2

Its estimator is

\sum_{i\in U}\sum_{j > i} \frac{(\pi_i\pi_j - \pi_{ij})}{\pi_{ij}} \Biggl(\frac{y_i}{\pi_i} - \frac{y_j}{\pi_j} \Biggr)^2

Examples

### Generate population data ---
N <- 500; n <- 50

set.seed(0)
x <- rgamma(500, scale=10, shape=5)
y <- abs( 2*x + 3.7*sqrt(x) * rnorm(N) )

pik  <- n * x/sum(x)
pikl <- jip_approx(pik, method='Hajek')

### Dummy sample ---
s   <- sample(N, n)


### Compute Variance ---
HTvar(y=y, pikl=pikl, sample=FALSE, method="HT")
HTvar(y=y, pikl=pikl, sample=FALSE, method="SYG")


### Estimate Variance ---
#' HTvar(y=y[s], pikl=pikl[s,s], sample=TRUE, method="HT")
#' HTvar(y=y[s], pikl=pikl[s,s], sample=TRUE, method="SYG")

rhobis/jipApprox documentation built on Sept. 12, 2023, 7:01 a.m.