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

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 (1999)Deville, 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 Deville (1999)Deville, 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..

A 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)

library(survey)

library(vardpoor)

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",

# poverty threshold range
order_quant = 50L ,

# 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
##  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 )
##  TRUE
# variances do not match exactly
attr( svypoormed( ~ eqincome , des_eusilc ) , 'var' )
##          eqincome
## eqincome  5311.47
varpoord_poormed_calculation$all_result$var
##  5302.086
# standard errors do not match exactly
varpoord_poormed_calculation$all_result$se
##  72.81542
SE( svypoormed( ~ eqincome , des_eusilc ) )
##          eqincome
## eqincome 72.87983

The variance estimate is computed by using the approximation defined in (1.1), where the linearized variable $$z$$ is defined by (1.2). 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
attr( svypoormed( ~ eqincome , des_eusilc_ultimate_cluster ) , 'var' )
##          eqincome
## eqincome 5302.086
varpoord_poormed_calculation$all_result$var
##  5302.086
# matches
varpoord_poormed_calculation$all_result$se
##  72.81542
SE( svypoormed( ~ eqincome , des_eusilc_ultimate_cluster ) )
##          eqincome
## eqincome 72.81542

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