For the most part, this document will present the functionalities of the function surveysd::calc.stError() which generates point estimates and standard errors for user-supplied estimation functions.

Prerequisites

In order to use a dataset with calc.stError(), several weight columns have to be present. Each weight column corresponds to a bootstrap sample. In the following examples, we will use the data from demo.eusilc() and attach the bootstrap weights using draw.bootstrap() and recalib(). Please refer to the documentation of those functions for more detail.

library(surveysd)

set.seed(1234)
eusilc <- demo.eusilc(prettyNames = TRUE)
dat_boot <- draw.bootstrap(eusilc, REP = 10, hid = "hid", weights = "pWeight",
                           strata = "region", period = "year")
dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region")
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps 
## 
## Convergence reached in 3 steps
##    year povertyRisk eqIncome onePerson  pWeight       w1       w2       w3
## 1: 2010       FALSE 16090.69     FALSE 504.5696 1.459066 1006.555 1.455580
## 2: 2010       FALSE 16090.69     FALSE 504.5696 1.459066 1006.555 1.455580
## 3: 2010       FALSE 16090.69     FALSE 504.5696 1.459066 1006.555 1.455580
## 4: 2011       FALSE 16090.69     FALSE 504.5696 1.380779  967.380 1.456154
## 5: 2011       FALSE 16090.69     FALSE 504.5696 1.380779  967.380 1.456154
##          w4       w5
## 1: 1.462730 1018.671
## 2: 1.462730 1018.671
## 3: 1.462730 1018.671
## 4: 1.516569 1004.087
## 5: 1.516569 1004.087

Estimator functions

The parameters fun and var in calc.stError() define the estimator to be used in the error analysis. There are two built-in estimator functions weightedSum() and weightedRatio() which can be used as follows.

povertyRate <- calc.stError(dat_boot_calib, var = "povertyRisk", fun = weightedRatio)
totalIncome <- calc.stError(dat_boot_calib, var = "eqIncome", fun = weightedSum)

Those functions calculate the ratio of persons at risk of povery (in percent) and the total income. By default, the results are calculated seperately for each reference period.

##    year     n       N val_povertyRisk stE_povertyRisk
## 1: 2010 14827 8182222        14.44422       0.5829182
## 2: 2011 14827 8182222        14.77393       0.6761383
## 3: 2012 14827 8182222        15.04515       0.4903267
## 4: 2013 14827 8182222        14.89013       0.4751962
## 5: 2014 14827 8182222        15.14556       0.5188642
## 6: 2015 14827 8182222        15.53640       0.4798736
## 7: 2016 14827 8182222        15.08315       0.3559527
## 8: 2017 14827 8182222        15.42019       0.5582783
##    year     n       N val_eqIncome stE_eqIncome
## 1: 2010 14827 8182222 162750998071    959510036
## 2: 2011 14827 8182222 161926931417    900201396
## 3: 2012 14827 8182222 162576509628   1195965007
## 4: 2013 14827 8182222 163199507862   1472304274
## 5: 2014 14827 8182222 163986275009   1525837944
## 6: 2015 14827 8182222 163416275447   1335217370
## 7: 2016 14827 8182222 162706205137   1204524680
## 8: 2017 14827 8182222 164314959107   1320123257

Columns that use the val_ prefix denote the point estimate belonging to the “main weight” of the dataset, which is pWeight in case of the dataset used here.

Columns with the stE_ prefix denote standard errors calculated with bootstrap replicates. The replicates result in using w1, w2, …, w10 instead of pWeight when applying the estimator.

n denotes the number of observations for the year and N denotes the total weight of those persons.

Custom estimators

In order to define a custom estimator function to be used in fun, the function needs to have two arguments like the example below.

## [1] TRUE

The parameters x and w can be assumed to be vectors with equal length with w being numeric and x being the column defined in the var argument. It will be called once for each period (in this case year) and for each weight column (in this case pWeight, w1, w2, …, w10).

Multiple estimators

In case an estimator should be applied to several columns of the dataset, var can be set to a vector containing all necessary columns.

multipleRates <- calc.stError(dat_boot_calib, var = c("povertyRisk", "onePerson"), fun = weightedRatio)
multipleRates$Estimates
##    year     n       N val_povertyRisk stE_povertyRisk val_onePerson
## 1: 2010 14827 8182222        14.44422       0.5829182      14.85737
## 2: 2011 14827 8182222        14.77393       0.6761383      14.85737
## 3: 2012 14827 8182222        15.04515       0.4903267      14.85737
## 4: 2013 14827 8182222        14.89013       0.4751962      14.85737
## 5: 2014 14827 8182222        15.14556       0.5188642      14.85737
## 6: 2015 14827 8182222        15.53640       0.4798736      14.85737
## 7: 2016 14827 8182222        15.08315       0.3559527      14.85737
## 8: 2017 14827 8182222        15.42019       0.5582783      14.85737
##    stE_onePerson
## 1:     0.2568455
## 2:     0.2774971
## 3:     0.2905952
## 4:     0.3452929
## 5:     0.4386546
## 6:     0.3803537
## 7:     0.3099259
## 8:     0.3025539

Here we see the relative number of persons at risk of poverty and the relative number of one-person households.

Grouping

The groups argument can be used to calculate estimators for different subsets of the data. This argument can take the grouping variable as a string that refers to a column name (usually a factor) in dat. If set, all estimators are not only split by the reference period but also by the grouping variable. For simplicity, only one reference period of the above data is used.

dat2 <- subset(dat_boot_calib, year == 2010)
for (att  in c("period", "weights", "b.rep"))
  attr(dat2, att) <- attr(dat_boot_calib, att)

To calculate the ratio of persons at risk of poverty for each federal state of austria, group = "region" can be used.

povertyRates <- calc.stError(dat2, var = "povertyRisk", fun = weightedRatio, group = "region")
povertyRates$Estimates
##     year     n       N        region val_povertyRisk stE_povertyRisk
##  1: 2010   549  260564    Burgenland        19.53984       1.7963347
##  2: 2010   733  377355    Vorarlberg        16.53731       3.2774567
##  3: 2010   924  535451      Salzburg        13.78734       2.2523362
##  4: 2010  1078  563648     Carinthia        13.08627       1.7427986
##  5: 2010  1317  701899         Tyrol        15.30819       1.3916058
##  6: 2010  2295 1167045        Styria        14.37464       1.3485703
##  7: 2010  2322 1598931        Vienna        17.23468       1.0419000
##  8: 2010  2804 1555709 Lower Austria        13.84362       1.7019212
##  9: 2010  2805 1421620 Upper Austria        10.88977       0.8701448
## 10: 2010 14827 8182222          <NA>        14.44422       0.5829182

The last column with region = NA denotes the aggregate over all regions. Note that the columns N and n now show the weighted and unweighted number of persons in each region.

Several grouping variables

In case more than one grouping variable is used, there are several options of calling calc.stError() depending on whether combinations of grouping levels should be regarded or not. We will consider the variables gender and region as our grouping variables and show three options on how calc.stError() can be called.

Option 1: All regions and all genders

Calculate the point estimate and standard error for each region and each gender. The number of rows in the output is therefore

\[n_\text{periods}\cdot(n_\text{regions} + n_\text{genders} + 1) = 1\cdot(9 + 2 + 1) = 12.\]

The last row is again the estimate for the whole period.

povertyRates <- calc.stError(dat2, var = "povertyRisk", fun = weightedRatio, 
                             group = c("gender", "region"))
povertyRates$Estimates
##     year     n       N gender        region val_povertyRisk
##  1: 2010   549  260564   <NA>    Burgenland        19.53984
##  2: 2010   733  377355   <NA>    Vorarlberg        16.53731
##  3: 2010   924  535451   <NA>      Salzburg        13.78734
##  4: 2010  1078  563648   <NA>     Carinthia        13.08627
##  5: 2010  1317  701899   <NA>         Tyrol        15.30819
##  6: 2010  2295 1167045   <NA>        Styria        14.37464
##  7: 2010  2322 1598931   <NA>        Vienna        17.23468
##  8: 2010  2804 1555709   <NA> Lower Austria        13.84362
##  9: 2010  2805 1421620   <NA> Upper Austria        10.88977
## 10: 2010  7267 3979572   male          <NA>        12.02660
## 11: 2010  7560 4202650 female          <NA>        16.73351
## 12: 2010 14827 8182222   <NA>          <NA>        14.44422
##     stE_povertyRisk
##  1:       1.7963347
##  2:       3.2774567
##  3:       2.2523362
##  4:       1.7427986
##  5:       1.3916058
##  6:       1.3485703
##  7:       1.0419000
##  8:       1.7019212
##  9:       0.8701448
## 10:       0.6294860
## 11:       0.6290441
## 12:       0.5829182

Option 2: All combinations of state and gender

Split the data by all cobinations of the two grouping variables. This will result in a larger output-table of the size

\[n_\text{periods}\cdot(n_\text{regions} \cdot n_\text{genders} + 1) = 1\cdot(9\cdot2 + 1)= 19.\]

povertyRates <- calc.stError(dat2, var = "povertyRisk", fun = weightedRatio, 
                             group = list(c("gender", "region")))
povertyRates$Estimates
##     year     n         N gender        region val_povertyRisk
##  1: 2010   261  122741.8   male    Burgenland       17.414524
##  2: 2010   288  137822.2 female    Burgenland       21.432598
##  3: 2010   359  182732.9   male    Vorarlberg       12.973259
##  4: 2010   374  194622.1 female    Vorarlberg       19.883637
##  5: 2010   440  253143.7   male      Salzburg        9.156964
##  6: 2010   484  282307.3 female      Salzburg       17.939382
##  7: 2010   517  268581.4   male     Carinthia       10.552149
##  8: 2010   561  295066.6 female     Carinthia       15.392924
##  9: 2010   650  339566.5   male         Tyrol       12.857542
## 10: 2010   667  362332.5 female         Tyrol       17.604861
## 11: 2010  1128  571011.7   male        Styria       11.671247
## 12: 2010  1132  774405.4   male        Vienna       15.590616
## 13: 2010  1167  596033.3 female        Styria       16.964539
## 14: 2010  1190  824525.6 female        Vienna       18.778813
## 15: 2010  1363  684272.5   male Upper Austria        9.074690
## 16: 2010  1387  772593.2 female Lower Austria       16.372949
## 17: 2010  1417  783115.8   male Lower Austria       11.348283
## 18: 2010  1442  737347.5 female Upper Austria       12.574206
## 19: 2010 14827 8182222.0   <NA>          <NA>       14.444218
##     stE_povertyRisk
##  1:       2.2233926
##  2:       2.0806293
##  3:       3.0976410
##  4:       3.7435376
##  5:       1.8974748
##  6:       2.5305822
##  7:       2.0566143
##  8:       1.9709443
##  9:       1.0350617
## 10:       2.2101768
## 11:       1.5095045
## 12:       1.3336364
## 13:       1.3743809
## 14:       0.9295877
## 15:       1.1700859
## 16:       1.8349672
## 17:       1.6423009
## 18:       0.7703255
## 19:       0.5829182

Option 3: Cobination of Option 1 and Option 2

In this case, the estimates and standard errors are calculated for

  • every gender,
  • every state and
  • every combination of state and gender.

The number of rows in the output is therefore

\[n_\text{periods}\cdot(n_\text{regions} \cdot n_\text{genders} + n_\text{regions} + n_\text{genders} + 1) = 1\cdot(9\cdot2 + 9 + 2 + 1) = 30.\]

povertyRates <- calc.stError(dat2, var = "povertyRisk", fun = weightedRatio, 
                             group = list("gender", "region", c("gender", "region")))
povertyRates$Estimates
##     year     n         N gender        region val_povertyRisk
##  1: 2010   261  122741.8   male    Burgenland       17.414524
##  2: 2010   288  137822.2 female    Burgenland       21.432598
##  3: 2010   359  182732.9   male    Vorarlberg       12.973259
##  4: 2010   374  194622.1 female    Vorarlberg       19.883637
##  5: 2010   440  253143.7   male      Salzburg        9.156964
##  6: 2010   484  282307.3 female      Salzburg       17.939382
##  7: 2010   517  268581.4   male     Carinthia       10.552149
##  8: 2010   549  260564.0   <NA>    Burgenland       19.539837
##  9: 2010   561  295066.6 female     Carinthia       15.392924
## 10: 2010   650  339566.5   male         Tyrol       12.857542
## 11: 2010   667  362332.5 female         Tyrol       17.604861
## 12: 2010   733  377355.0   <NA>    Vorarlberg       16.537310
## 13: 2010   924  535451.0   <NA>      Salzburg       13.787343
## 14: 2010  1078  563648.0   <NA>     Carinthia       13.086268
## 15: 2010  1128  571011.7   male        Styria       11.671247
## 16: 2010  1132  774405.4   male        Vienna       15.590616
## 17: 2010  1167  596033.3 female        Styria       16.964539
## 18: 2010  1190  824525.6 female        Vienna       18.778813
## 19: 2010  1317  701899.0   <NA>         Tyrol       15.308190
## 20: 2010  1363  684272.5   male Upper Austria        9.074690
## 21: 2010  1387  772593.2 female Lower Austria       16.372949
## 22: 2010  1417  783115.8   male Lower Austria       11.348283
## 23: 2010  1442  737347.5 female Upper Austria       12.574206
## 24: 2010  2295 1167045.0   <NA>        Styria       14.374637
## 25: 2010  2322 1598931.0   <NA>        Vienna       17.234683
## 26: 2010  2804 1555709.0   <NA> Lower Austria       13.843623
## 27: 2010  2805 1421620.0   <NA> Upper Austria       10.889773
## 28: 2010  7267 3979571.7   male          <NA>       12.026600
## 29: 2010  7560 4202650.3 female          <NA>       16.733508
## 30: 2010 14827 8182222.0   <NA>          <NA>       14.444218
##     year     n         N gender        region val_povertyRisk
##     stE_povertyRisk
##  1:       2.2233926
##  2:       2.0806293
##  3:       3.0976410
##  4:       3.7435376
##  5:       1.8974748
##  6:       2.5305822
##  7:       2.0566143
##  8:       1.7963347
##  9:       1.9709443
## 10:       1.0350617
## 11:       2.2101768
## 12:       3.2774567
## 13:       2.2523362
## 14:       1.7427986
## 15:       1.5095045
## 16:       1.3336364
## 17:       1.3743809
## 18:       0.9295877
## 19:       1.3916058
## 20:       1.1700859
## 21:       1.8349672
## 22:       1.6423009
## 23:       0.7703255
## 24:       1.3485703
## 25:       1.0419000
## 26:       1.7019212
## 27:       0.8701448
## 28:       0.6294860
## 29:       0.6290441
## 30:       0.5829182
##     stE_povertyRisk