Stratified survey calculations by hand and with survey package don't agree. Simulation results

Question

Bounty info: I originally emailed Thomas Lumley at an old email address. He did respond to an email to his new address.

Note: Long post (lots of code)

I can’t seem to replicate the results of the survey function using very basic by-hand calculations. I’m also having trouble understanding some aspects of the survey function.

I created a simulation to test this in a controlled way.

The formulas I used for the by-hand calculations are at the end.

Create Population

    set.seed(05022020)
    ### Total size
    P = 1000

    ### mean
    mu = 10

    ### sd
    sigma = 5

    population_data = rnorm(n=P, mean=mu, sd=sigma)

Stratify the real population

I want to pretend we have 5 different groups in the population that are more homogeneous inside them than between them.

The algorithm I will use is

1. Sort the values
2. Randomly pick 4 numbers between 0 and P
3. Use these as cutoffs for the strata

    cutoffs = sample(0:P, 4) %>% sort()
    print(cutoffs)
    #> [1]  20 156 564 868

    population =

      ### Start with data
      data.frame(
        data = sort(population_data),
        ### Unit ID
        UID = paste0("Unit_", 1:P)
        ) %>%

      ### Create strata
      mutate(stratum=
        case_when(
          row_number() <= cutoffs[1] ~ "Stratum_1",
          row_number() > cutoffs[1] & row_number() <= cutoffs[2] ~ "Stratum_2",
          row_number() > cutoffs[2] & row_number() <= cutoffs[3] ~ "Stratum_3",
          row_number() > cutoffs[3] & row_number() <=cutoffs[4] ~ "Stratum_4",
          TRUE ~ "Stratum_5"
        ) %>%
        factor(levels=paste0('Stratum_', 1:5))
      ) %>%

      ### calculate strata populations
      group_by(stratum) %>%
      mutate(stratum_population = n()) %>%
      ungroup()

    summary(population)
    #>       data                UID           stratum    stratum_population
    #>  Min.   :-12.658   Unit_1   :  1   Stratum_1: 20   Min.   : 20.0     
    #>  1st Qu.:  6.424   Unit_10  :  1   Stratum_2:136   1st Qu.:136.0     
    #>  Median : 10.119   Unit_100 :  1   Stratum_3:408   Median :304.0     
    #>  Mean   :  9.938   Unit_1000:  1   Stratum_4:304   Mean   :295.2     
    #>  3rd Qu.: 13.499   Unit_101 :  1   Stratum_5:132   3rd Qu.:408.0     
    #>  Max.   : 25.062   Unit_102 :  1                   Max.   :408.0     
    #>                    (Other)  :994

    # population %>%
    #   ggplot(aes(x=stratum, y=data, color=stratum)) +
    #   geom_beeswarm(alpha=0.3, size=1, stroke=F) +
    #   geom_boxplot(alpha=0)

Note: The actual sample mean of the entire population is 9.94.

Stratified Survey

Pick SRSWOR

Note: Because samples are chosen at random, this is approximately population-weighted (self-weighted).

    # set.seed(05022020)

    sample_size = 100

    srswor =
      population %>%
      sample_n(sample_size)

    srswor %>%
      ggplot(aes(x=stratum, y=data, color=stratum)) +
      geom_beeswarm(alpha=0.3, size=1, stroke=F) +
      geom_boxplot(alpha=0)

Raw Stats

    srswor %>%
      summarize(
        mu=mean(data),
        sd=sd(data),
        s2=var(data),
        SE=sqrt(s2/n()),
        CI_low = mu - qnorm(0.975)*SE,
        CI_high = mu + qnorm(0.975)*SE
      )
    #> # A tibble: 1 x 6
    #>      mu    sd    s2    SE CI_low CI_high
    #>   <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>
    #> 1  9.54  5.27  27.8 0.527   8.51    10.6

Unstratified calculation with FPC

By Hand

    unstratified_summary =
      srswor %>%
      ### Collect global statistics
      mutate(
        ### N: total population
        N = P,
        ### H: number of strata
        H = length(unique(stratum)),
        ### n: total number of samples
        n = n(), # number of rows
        ### \hat mu: population mean
        mu = mean(data),
        ### \hat s
        s2 = var(data),
        ### FPC: finite population correction
        ###      entire population for SRSWOR
        FPC = (N-n)/N
        # FPC = 1
      ) %>%
      ### Just keep needed columns
      select(N, n, H, mu, s2, FPC) %>%
      ### Just keep one unique row for the entire population
      unique()

    print(unstratified_summary)
    #> # A tibble: 1 x 6
    #>       N     n     H    mu    s2   FPC
    #>   <dbl> <int> <int> <dbl> <dbl> <dbl>
    #> 1  1000   100     5  9.54  27.8   0.9

    unstratified_stats =
      unstratified_summary %>%
      mutate(
        sd = sqrt(s2),
        SE = sqrt(s2/n*FPC),
        CI_low = mu - qnorm(0.975)*SE,
        CI_high = mu + qnorm(0.975)*SE
      ) %>%
      select(mu, sd, SE, CI_low, CI_high)

    print(unstratified_stats)
    #> # A tibble: 1 x 5
    #>      mu    sd    SE CI_low CI_high
    #>   <dbl> <dbl> <dbl>  <dbl>   <dbl>
    #> 1  9.54  5.27 0.500   8.56    10.5

Note: Using the FPC did narrow the CI as expected.

Using Survey Function

    unstrat_design = svydesign(
      id = ~ 1,
      strata=NULL,
      FPC = ~FPC,
      data = srswor %>% mutate(FPC = 1 - sample_size/P)
      # data = srswor
    )
    #> Warning in svydesign.default(id = ~1, strata = NULL, FPC = ~FPC, data = srswor
    #> %>% : No weights or probabilities supplied, assuming equal probability
    print('')
    #> [1] ""
    unstrat_mean = svymean(~data, unstrat_design)
    print(unstrat_mean)
    #>        mean     SE
    #> data 9.5438 0.5273
    confint(unstrat_mean)
    #>         2.5 %   97.5 %
    #> data 8.510296 10.57725

Note:

1. It looks like svymean is calling the $\hat \sigma$ the SE.
2. The svymean function gives a wider CI.
3. svymean gave the same calculation as raw, i.e. without using the FPC provided.

Stratified Calculations

By Hand

    stratified_summary =
      srswor %>%
      ### Collect global statistics
      mutate(
        ### N: POPULATION size
        N = P, # defined size of population
        ### H: number of strata
        H = length(unique(stratum)),
        ### n: total number of samples, sample size
        n = n(), # number of rows in sample
        ### \hat mu: population mean
        mu = mean(data),
        ### \hat s
        s2 = var(data),
        ### FPC: finite population correction
        ###      entire population for SRSWOR
        FPC = (N-n)/N
      ) %>%

      ### Collect strata statistics
      group_by(stratum) %>%
      mutate(
        ### N_h: POPULATION stratum size
        ###      calculated when srswor s defined aboev.
        N_h = stratum_population, 
        ### \hat mu_h: stratum sample mean
        mu_h = mean(data),
        ### n_h: SAMPLE stratum size
        n_h = n(),
        ### \hat s^2_h: stratum sample variance
        s2_h = var(data),
        ### Weight
        w_h = N_h/N, # sum w_h = 1
        ### fpc: strata level for SRSWOR
        # fpc = (N - N_h)/(N-1)
        fpc = 1 - w_h
      ) %>%
      ungroup() %>%

      ### Just keep things we need
      select(stratum, N_h, n_h, mu_h, s2_h, w_h, fpc, FPC) %>%

      ### Just keep one unique row per stratum
      unique()

    print(sprintf("number of rows: %d", nrow(stratified_summary)))
    #> [1] "number of rows: 5"
    head(stratified_summary)
    #> # A tibble: 5 x 8
    #>   stratum     N_h   n_h  mu_h  s2_h   w_h   fpc   FPC
    #>   <fct>     <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
    #> 1 Stratum_4   304    22 13.2   2.13 0.304 0.696   0.9
    #> 2 Stratum_2   136    19  2.40  2.01 0.136 0.864   0.9
    #> 3 Stratum_3   408    44  8.32  3.26 0.408 0.592   0.9
    #> 4 Stratum_5   132    14 18.1   4.28 0.132 0.868   0.9
    #> 5 Stratum_1    20     1 -1.21 NA    0.02  0.98    0.9

    stratified_stats =
      stratified_summary %>%
      drop_na() %>%
      filter(n_h>1) %>%
      summarize(
        hat_mu = sum(w_h*mu_h),
        ### Note: `na.rm=T` takes care of the situation where
        ### n_h = 1
        hat_se = sqrt(sum(w_h^2*s2_h/n_h*fpc)),
        CI_low = hat_mu - qnorm(0.975)*hat_se,
        CI_high = hat_mu + qnorm(0.975)*hat_se
      )

    print(stratified_stats)
    #> # A tibble: 1 x 4
    #>   hat_mu hat_se CI_low CI_high
    #>    <dbl>  <dbl>  <dbl>   <dbl>
    #> 1   10.1  0.141   9.85    10.4

Using Survey Function

Combine srswor data with summary stats to

1. Filter for $n\_h > 1$
2. Have a column with fpc

    temp_data =
      inner_join(stratified_summary, srswor, by="stratum") %>%
      select(stratum, data, stratum_population, n_h, fpc)

    dim(temp_data)
    #> [1] 100   5

    summary(temp_data)
    #>       stratum        data        stratum_population      n_h       
    #>  Stratum_1: 1   Min.   :-1.208   Min.   : 20.0      Min.   : 1.00  
    #>  Stratum_2:19   1st Qu.: 5.447   1st Qu.:136.0      1st Qu.:19.00  
    #>  Stratum_3:44   Median : 9.441   Median :304.0      Median :22.00  
    #>  Stratum_4:22   Mean   : 9.544   Mean   :290.9      Mean   :29.78  
    #>  Stratum_5:14   3rd Qu.:13.024   3rd Qu.:408.0      3rd Qu.:44.00  
    #>                 Max.   :24.229   Max.   :408.0      Max.   :44.00  
    #>       fpc        
    #>  Min.   :0.5920  
    #>  1st Qu.:0.5920  
    #>  Median :0.6960  
    #>  Mean   :0.7091  
    #>  3rd Qu.:0.8640  
    #>  Max.   :0.9800

    head(temp_data)
    #> # A tibble: 6 x 5
    #>   stratum    data stratum_population   n_h   fpc
    #>   <fct>     <dbl>              <int> <int> <dbl>
    #> 1 Stratum_4  10.9                304    22 0.696
    #> 2 Stratum_4  14.6                304    22 0.696
    #> 3 Stratum_4  14.1                304    22 0.696
    #> 4 Stratum_4  14.9                304    22 0.696
    #> 5 Stratum_4  12.1                304    22 0.696
    #> 6 Stratum_4  14.6                304    22 0.696

Attempt 1

Use the stratum population for the FPC as per the documentation

    data_design = svydesign(
      id = ~ 1,
      fpc = ~ stratum_population,
      strata = ~ stratum,
      data = temp_data %>% filter(n_h> 1)
    )
    print('')
    #> [1] ""

    svymean(~data, data_design)
    #>        mean     SE
    #> data 10.329 0.1637

    confint(svymean(~data, data_design))
    #>         2.5 %  97.5 %
    #> data 10.00855 10.6502

Note:
1. My mean is closer to the real one.
2. My CI is narrower.

Attempt 2

Add in the FPC as calculated above

    data_design = svydesign(
      id = ~ 1,
      fpc = ~ fpc,
      strata = ~ stratum,
      data = temp_data %>% filter(n_h> 1)
    )
    print('')
    #> [1] ""

    svymean(~data, data_design)
    #>        mean     SE
    #> data 9.5823 0.1015

    confint(svymean(~data, data_design))
    #>         2.5 %   97.5 %
    #> data 9.383343 9.781196

Note: I don't get this at all.

^{Created on 2020-05-05 by the reprex package (v0.3.0)}

Formulas

True Known Values

$H$ --- number of strata

$N_h$ --- number of units in stratum $h$. (not always known) $\sum_{h \in \mathcal{S}} N_h = N_\mathcal{S}$

$n_h$ --- number of samples, i.e. units actually sampled in stratum $h$. $\sum_{h \in \mathcal{S}} n_h = n_\mathcal{S}$

$\mu_h = \overline y_h$ --- unknown true mean of $y$ in stratum $h$

$\tau_h = \sum_{y \in h} y = N_h \mu_h$ --- unknown true total of $y$ in stratum $h$

$\sum_{h \in \mathcal{S}} \tau_h = \sum_{h \in \mathcal{S}} N_h \mu_h = \tau_\mathcal{S}$ --- total of $y$ in population

$\sigma_h^2 = \overline {(y - \overline y_h)^2}$ --- unknown true variance of $y$ in stratum $h$

Estimated Values

$\hat \mu_h = \widehat {\overline {y_h}}$ --- sample mean of stratum $h$, unbiased estimator of $\mu_h$

$\hat \sigma_h^2 = s_h^2 = \frac{1}{n_h-1}\sum_{\hat h} (y - \widehat {\overline {y_h}})^2 = \widehat{\overline{(y - \widehat {\overline {y_h}})^2}}$ --- sample variance of stratum $h$, unbiased estimate of $\sigma_h$

$\hat \tau_h = N_h \hat \mu_h$ --- unbiased estimate of $\tau_h$

$\widehat{\mathrm{V}}(\hat \mu_h) = \frac{\hat{\sigma}^2}{n_h}$ --- unbiased estimator of the variance of the stratum mean for SRSWR

$\widehat{\mathrm{V}}(\hat \mu_h) = \frac{\hat{\sigma}^2}{n_h}\frac{N_h - n_h}{N_h}$ --- unbiased estimator of the variance of the stratum mean for SRSWOR

$\frac{N_h - n_h}{N_h}$ --- adjustment factor for sampling without replacement (SRSWOR)

$w_h = \frac{n_h}{N_h}$ --- stratum weight (self weight)

Population Estimates from Strata Statistics

Unbiased estimator for the population total

$$\mathrm{E}(\hat \tau_\mathcal{S}) = \sum_{h \in \mathcal{S}} \hat \tau_h = \sum_{h \in \mathcal{S}} N_h \hat \mu_h$$

Unbiased estimator for the population mean

$$\mathrm{E}(\hat \mu) = \sum_{h \in \mathcal{S}} w_h\hat \mu_h$$

Unbiased estimator for the population variance

$$\mathrm{E}(\hat \sigma^2) = \sum_{h \in \mathcal{S}} \hat \sigma_h^2$$

Unbiased estimator for the variance of the population mean (SRSWOR)

$$\widehat{\mathrm{V}}(\hat \mu) = \sum_{h \in \mathcal{S}} w_h^2 \frac{\hat \sigma^2_h}{n_h}\frac{N_h - n_h}{N_h}$$

Unbiased estimator for the variance of the population total (SRSWOR)

$$\widehat{\mathrm{V}}(\hat \tau) = \sum_{h \in \mathcal{S}} N_h^2 \frac{\hat \sigma^2_h}{n_h}\frac{N_h - n_h}{N_h}$$

Answer 1

Let's start with this

unstrat_design = svydesign(
      id = ~ 1,
      strata=NULL,
      FPC = ~FPC,
      data = srswor %>% mutate(FPC = 1-sample_size/P)
    )
## Warning in svydesign.default(id = ~1, strata = NULL, FPC = ~FPC, data
## = srswor %>% : No weights or probabilities supplied, assuming equal
## probability
unstrat_design
## Independent Sampling design (with replacement)
## svydesign(id = ~1, strata = NULL, FPC = ~FPC, data = srswor %>% 
##     mutate(FPC = 1 - sample_size/P))
correct_unstrat_design = svydesign(
      id = ~ 1,
      strata=NULL,
      fpc = ~FPC,
      data = srswor %>% mutate(FPC = sample_size/P)
    )

correct_unstrat_design
## Independent Sampling design
## svydesign(id = ~1, strata = NULL, fpc = ~FPC, data = srswor %>% 
##     mutate(FPC = sample_size/P))
unstrat_mean = svymean(~data, correct_unstrat_design)
print(unstrat_mean)
##        mean     SE
## data 9.5438 0.5002

svydesign does not have an FPC argument; it has an fpc argument. Unfortunately, the way the S3 method/inheritance system is set up, a method must accept and ignore arguments that it does not recognise.

Also, the help page for svydesign says

The finite population correction can be specified either as the total population size in each stratum or as the fraction of the total population that has been sampled. In either case the relevant population size is the sampling units. That is, sampling 100 units from a population stratum of size 500 can be specified as 500 or as 100/500=0.2. The exception is for PPS sampling without replacement, where the sampling probability (which will be different for each PSU) must be used.

That is, you can use fpc=sample_size or fpc=sample_size/P, but not fpc=1-sample_size/P

There are some situations where svymean does give different results from what a textbook might suggest. These are situations where the population size is known, but svydesign has not been told that it is known. The most obvious example is cluster sampling, where a user might know the sizes of unsampled clusters and thus the population size, and a textbook might estimate the mean as the estimated total divided by the known population size, but svydesign will divide the estimated total by the estimated population size. In that situation, calibrate() can be used to provide the extra information.

Even when svymean gives the same answer, it doesn't necessarily use the same formula, because it has to work more generally.

I'll go on to look at the stratified example, but this gets us a lot of the way.

Actually, the stratified sample has got a bit confusing -- for a start, it isn't actually a stratified sample. What I will do is analyse a stratified sample that we already have available in the package

library(foreign)
write.dta(model.frame(correct_unstrat_design),file="mean.dta")
data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
dstrat
## Stratified Independent Sampling design
## svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
##     fpc = ~fpc)
svymean(~enroll, dstrat)
##          mean     SE
## enroll 595.28 18.509

Now I'll do the same calculation a couple of different ways in R: the mean is also a ratio estimator and is also a regression coefficient

summary(svyglm(enroll~1, dstrat))
## 
## Call:
## svyglm(formula = enroll ~ 1, design = dstrat)
## 
## Survey design:
## svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
##     fpc = ~fpc)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   595.28      18.51   32.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 195567.7)
## 
## Number of Fisher Scoring iterations: 2
dstrat<-update(dstrat, one=1+0*enroll)
svyratio(~enroll, ~one, dstrat)
## Ratio estimator: svyratio.survey.design2(~enroll, ~one, dstrat)
## Ratios=
##             one
## enroll 595.2821
## SEs=
##             one
## enroll 18.50851

And I'll go and do it in Stata and get the same answer

. use apistrat

. svyset snum [pw=pw], fpc(fpc) str(stype)

      pweight: pw
          VCE: linearized
  Single unit: missing
     Strata 1: stype
         SU 1: snum
        FPC 1: fpc

. svy: mean enroll
(running mean on estimation sample)

Survey: Mean estimation

Number of strata =       3        Number of obs   =        200
Number of PSUs   =     200        Population size =      6,194
                                  Design df       =        197

--------------------------------------------------------------
             |             Linearized
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
      enroll |   595.2821   18.50851      558.7819    631.7824
--------------------------------------------------------------

Next, does this agree with the hand calculations?

(muhat_h<-with(apistrat, by(enroll, stype, mean)))
## stype: E
## [1] 416.78
## -------------------------------------------------------- 
## stype: H
## [1] 1320.7
## -------------------------------------------------------- 
## stype: M
## [1] 832.48
(s2hat_h<-with(apistrat, by(enroll, stype, var)))
## stype: E
## [1] 27576.88
## -------------------------------------------------------- 
## stype: H
## [1] 450339.9
## -------------------------------------------------------- 
## stype: M
## [1] 156307.3
(Nh <-with(apipop, by(enroll, stype,length)))
## stype: E
## [1] 4421
## -------------------------------------------------------- 
## stype: H
## [1] 755
## -------------------------------------------------------- 
## stype: M
## [1] 1018
(nh<-with(apistrat, by(enroll, stype,length)))
## stype: E
## [1] 100
## -------------------------------------------------------- 
## stype: H
## [1] 50
## -------------------------------------------------------- 
## stype: M
## [1] 50
(tauhat_h <- Nh*muhat_h)
## stype: E
## [1] 1842584
## -------------------------------------------------------- 
## stype: H
## [1] 997128.5
## -------------------------------------------------------- 
## stype: M
## [1] 847464.6
(Vhatmuhat_h <- s2hat_h/nh)
## stype: E
## [1] 275.7688
## -------------------------------------------------------- 
## stype: H
## [1] 9006.799
## -------------------------------------------------------- 
## stype: M
## [1] 3126.146
(Vhatmuhat_again_h <- (s2hat_h/nh)*(Nh-nh)/nh)
## stype: E
## [1] 11915.97
## -------------------------------------------------------- 
## stype: H
## [1] 126995.9
## -------------------------------------------------------- 
## stype: M
## [1] 60522.18
(w_h <- nh/Nh)
## stype: E
## [1] 0.02261932
## -------------------------------------------------------- 
## stype: H
## [1] 0.06622517
## -------------------------------------------------------- 
## stype: M
## [1] 0.04911591
(Emuhat=sum(w_h*muhat_h))
## [1] 137.7789
(Vhatmuhat_yetagain_h = sum(w_h^2*(s2hat_h/nh)*(Nh-nh)/Nh))
## [1] 44.19468
sqrt(Vhatmuhat_yetagain_h)
## [1] 6.647908

Apparently not. One problem is the definition of $w_h=n_h/N_h$. This has to be wrong in the equation for the unbiased estimate of the mean, because it would make the mean smaller if the sampling fraction were smaller. In that equation it should be $w_h=N_h/\sum_h N_h$: the population fraction that each stratum makes up.

With that correction

w_h<- Nh/sum(Nh)
(Emuhat=sum(w_h*muhat_h))
## [1] 595.2821
(Vhatmuhat_yetagain_h = sum(w_h^2*(s2hat_h/nh)*(Nh-nh)/Nh))
## [1] 342.565
sqrt(Vhatmuhat_yetagain_h)
## [1] 18.50851

And it matches!

And, last of all, the Maintainer email on a CRAN package is required to be up to date and to deliver to the package maintainer. Some packages have a separate issue tracker, and you should obviously use that if it's given; others (in particular, those older than github) may not.