3.5 Median Income Below the At Risk of Poverty Threshold (svypoormed)

✔️ median income among those below the threshold
✔️ useful for understanding the depth of poverty
✔️ related to the RMPG
❌ not very common outside the EU
❌ not immediately interpretable in terms of the poverty gap

Median income below the at-risk-of-poverty-threshold (POORMED) is median of incomes of people having the income below the ARPT:

\[ poormed = median\{y_i; y_i< arpt\} \] The details of the linearization of the POORMED are discussed by Deville (⊕1999Deville, Jean-Claude. 1999. “Variance Estimation for Complex Statistics and Estimators: Linearization and Residual Techniques.” Survey Methodology 25 (2): 193–203. http://www.statcan.gc.ca/pub/12-001-x/1999002/article/4882-eng.pdf.) and Osier (⊕2009Osier, Guillaume. 2009. “Variance Estimation for Complex Indicators of Poverty and Inequality.” Journal of the European Survey Research Association 3 (3): 167–95. http://ojs.ub.uni-konstanz.de/srm/article/view/369.).

3.5.1 Replication Example

The R vardpoor package (⊕Breidaks, Liberts, and Ivanova 2016Breidaks, Juris, Martins Liberts, and Santa Ivanova. 2016. “Vardpoor: Estimation of Indicators on Social Exclusion and Poverty and Its Linearization, Variance Estimation.” Riga, Latvia: CSB.), created by researchers at the Central Statistical Bureau of Latvia, includes a POORMED coefficient calculation using the ultimate cluster method. The example below reproduces those statistics.

Load and prepare the same data set:

# load the convey package
library(convey)

# load the survey library
library(survey)

# load the vardpoor library
library(vardpoor)

# load the vardpoor library
library(laeken)

# load the synthetic EU statistics on income & living conditions
data(eusilc)

# make all column names lowercase
names(eusilc) <- tolower(names(eusilc))

# add a column with the row number
dati <- data.table::data.table(IDd = 1:nrow(eusilc), eusilc)

# calculate the poormed coefficient
# using the R vardpoor library
varpoord_poormed_calculation <-
  varpoord(
    # analysis variable
    Y = "eqincome",
    
    # weights variable
    w_final = "rb050",
    
    # row number variable
    ID_level1 = "IDd",
    
    # row number variable
    ID_level2 = "IDd",
    
    # strata variable
    H = "db040",
    
    N_h = NULL ,
    
    # clustering variable
    PSU = "rb030",
    
    # data.table
    dataset = dati,
    
    # poormed coefficient function
    type = "linpoormed",
    
    # get linearized variable
    outp_lin = TRUE
    
  )



# construct a survey.design
# using our recommended setup
des_eusilc <-
  svydesign(
    ids = ~ rb030 ,
    strata = ~ db040 ,
    weights = ~ rb050 ,
    data = eusilc
  )

# immediately run the convey_prep function on it
des_eusilc <- convey_prep(des_eusilc)

# coefficients do match
varpoord_poormed_calculation$all_result$value

## [1] 8803.735

coef(svypoormed( ~ eqincome , des_eusilc))

## eqincome 
## 8803.735

# linearized variables do match
# vardpoor
lin_poormed_varpoord <-
  varpoord_poormed_calculation$lin_out$lin_poormed
# convey
lin_poormed_convey <-
  attr(svypoormed( ~ eqincome , des_eusilc), "lin")

# check equality
all.equal(lin_poormed_varpoord, lin_poormed_convey)

## [1] "names for current but not for target"

# variances do not match exactly
attr(svypoormed( ~ eqincome , des_eusilc) , 'var')

##          eqincome
## eqincome  5311.47

varpoord_poormed_calculation$all_result$var

## [1] 5302.086

# standard errors do not match exactly
varpoord_poormed_calculation$all_result$se

## [1] 72.81542

SE(svypoormed( ~ eqincome , des_eusilc))

##          eqincome
## eqincome 72.87983

The variance estimator and the linearized variable \(z\) are both defined in Linearization-Based Variance Estimation. The functions convey::svypoormed and vardpoor::linpoormed produce the same linearized variable \(z\).

However, the measures of uncertainty do not line up, because library(vardpoor) defaults to an ultimate cluster method that can be replicated with an alternative setup of the survey.design object.

# within each strata, sum up the weights
cluster_sums <-
  aggregate(eusilc$rb050 , list(eusilc$db040) , sum)

# name the within-strata sums of weights the `cluster_sum`
names(cluster_sums) <- c("db040" , "cluster_sum")

# merge this column back onto the data.frame
eusilc <- merge(eusilc , cluster_sums)

# construct a survey.design
# with the fpc using the cluster sum
des_eusilc_ultimate_cluster <-
  svydesign(
    ids = ~ rb030 ,
    strata = ~ db040 ,
    weights = ~ rb050 ,
    data = eusilc ,
    fpc = ~ cluster_sum
  )

# again, immediately run the convey_prep function on the `survey.design`
des_eusilc_ultimate_cluster <-
  convey_prep(des_eusilc_ultimate_cluster)


# matches
stopifnot(all.equal(
  attr(svypoormed(~ eqincome , des_eusilc_ultimate_cluster) , 'var')[1],
  varpoord_poormed_calculation$all_result$var
))

# matches
stopifnot(all.equal(
  SE(svypoormed(~ eqincome , des_eusilc_ultimate_cluster))[1],
  varpoord_poormed_calculation$all_result$se
))

For additional usage examples of svypoormed, type ?convey::svypoormed in the R console.

3.5.2 Real World Examples

This section displays example results using nationally-representative surveys from both the United States and Brazil. We present a variety of surveys, levels of analysis, and subpopulation breakouts to provide users with points of reference for the range of plausible values of the svypoormed function.

To understand the construction of each survey design object and respective variables of interest, please refer to section 1.4 for CPS-ASEC, section 1.5 for PNAD Contínua, and section 1.6 for SCF.

3.5.2.1 CPS-ASEC Household Income

svypoormed( ~ htotval , cps_household_design)

##         poormed     SE
## htotval   24001 202.67

svyby( ~ htotval , ~ sex , cps_household_design , svypoormed)

##           sex htotval se.htotval
## male     male   25010   248.9748
## female female   23040   376.7984

3.5.2.2 CPS-ASEC Family Income

svypoormed(~ ftotval , cps_family_design)

##         poormed     SE
## ftotval   29000 221.99

svyby(~ ftotval , ~ sex , cps_family_design , svypoormed)

##           sex ftotval se.ftotval
## male     male   30044   275.2839
## female female   27184   496.5440

3.5.2.3 CPS-ASEC Worker Earnings

svypoormed( ~ pearnval , cps_ftfy_worker_design)

##          poormed     SE
## pearnval   28000 684.82

svyby( ~ pearnval , ~ sex , cps_ftfy_worker_design , svypoormed)

##           sex pearnval se.pearnval
## male     male    28000    144.2578
## female female    28000    878.5112

3.5.2.4 PNAD Contínua Per Capita Income

svypoormed( ~ deflated_per_capita_income , pnadc_design , na.rm = TRUE)

##                            poormed     SE
## deflated_per_capita_income  366.23 4.4574

svyby( ~ deflated_per_capita_income ,
       ~ sex ,
       pnadc_design ,
       svypoormed ,
       na.rm = TRUE)

##           sex deflated_per_capita_income se.deflated_per_capita_income
## male     male                   369.2743                      4.504995
## female female                   364.3737                      4.483770

3.5.2.5 PNAD Contínua Worker Earnings

svypoormed( ~ deflated_labor_income , pnadc_design , na.rm = TRUE)

##                       poormed     SE
## deflated_labor_income  502.67 4.5457

svyby( ~ deflated_labor_income , ~ sex , pnadc_design , svypoormed , na.rm = TRUE)

##           sex deflated_labor_income se.deflated_labor_income
## male     male              512.4121                5.2748713
## female female              498.4078                0.8745423

3.5.2.6 SCF Family Net Worth

scf_MIcombine(with(scf_design , svypoormed( ~ networth)))

## Multiple imputation results:
##       m <- length(results)
##       scf_MIcombine(with(scf_design, svypoormed(~networth)))
##          results       se
## networth   14242 963.1121

scf_MIcombine(with(scf_design , svyby( ~ networth, ~ hhsex , svypoormed)))

## Multiple imputation results:
##       m <- length(results)
##       scf_MIcombine(with(scf_design, svyby(~networth, ~hhsex, svypoormed)))
##        results       se
## male   19082.0 1907.857
## female  8896.4 1268.692

3.5.2.7 SCF Family Income

scf_MIcombine(with(scf_design , svypoormed( ~ income)))

## Multiple imputation results:
##       m <- length(results)
##       scf_MIcombine(with(scf_design, svypoormed(~income)))
##         results       se
## income 25077.15 887.2537

scf_MIcombine(with(scf_design , svyby( ~ income, ~ hhsex , svypoormed)))

## Multiple imputation results:
##       m <- length(results)
##       scf_MIcombine(with(scf_design, svyby(~income, ~hhsex, svypoormed)))
##        results       se
## male   27022.8 702.1986
## female 23347.7 961.0849