The R code below replicates multiple statistics and standard errors from the presentation by OPHI Research Officer Christop Jindra:
library(survey)
library(convey)
library(webuse)
# load the same microdata set used by Jindra in his presentation
webuse("nlsw88")
# coerce that `tbl_df` to a standard R `data.frame`
nlsw88 <- data.frame(nlsw88)
# create a `collgrad` column
nlsw88$collgrad <-
factor(
as.numeric(nlsw88$collgrad) ,
label = c('not college grad' , 'college grad') ,
ordered = TRUE
)
# initiate a linearization-based survey design object
des_nlsw88 <- svydesign(ids = ~ 1 , data = nlsw88)## Warning in svydesign.default(ids = ~1, data = nlsw88): No weights or
## probabilities supplied, assuming equal probability# immediately run the `convey_prep` function on the survey design
des_nlsw88 <- convey_prep(des_nlsw88)
# PDF page 9
result <-
svyafc(
~ wage + collgrad + hours,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = 1 / 3,
g = 0,
na.rm = T
)
result## alkire-foster SE
## [1,] 0.36991 0.0053## coef SE
## H 0.8082070 0.008316807
## A 0.4576895 0.004573443# PDF page 10
for (ks in seq(0.1, 1, .1)) {
result <-
svyafc(
~ wage + collgrad + hours,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = ks ,
g = 0,
na.rm = T
)
print(result)
print(attr(result , 'extra'))
}## alkire-foster SE
## [1,] 0.36991 0.0053
## coef SE
## H 0.8082070 0.008316807
## A 0.4576895 0.004573443
## alkire-foster SE
## [1,] 0.36991 0.0053
## coef SE
## H 0.8082070 0.008316807
## A 0.4576895 0.004573443
## alkire-foster SE
## [1,] 0.36991 0.0053
## coef SE
## H 0.8082070 0.008316807
## A 0.4576895 0.004573443
## alkire-foster SE
## [1,] 0.18659 0.0068
## coef SE
## H 0.2582516 0.009245464
## A 0.7225101 0.005174483
## alkire-foster SE
## [1,] 0.18659 0.0068
## coef SE
## H 0.2582516 0.009245464
## A 0.7225101 0.005174483
## alkire-foster SE
## [1,] 0.18659 0.0068
## coef SE
## H 0.2582516 0.009245464
## A 0.7225101 0.005174483
## alkire-foster SE
## [1,] 0.043265 0.0043
## coef SE
## H 0.04326494 0.004297766
## A 1.00000000 0.000000000
## alkire-foster SE
## [1,] 0.043265 0.0043
## coef SE
## H 0.04326494 0.004297766
## A 1.00000000 0.000000000
## alkire-foster SE
## [1,] 0.043265 0.0043
## coef SE
## H 0.04326494 0.004297766
## A 1.00000000 0.000000000
## alkire-foster SE
## [1,] 0.043265 0.0043
## coef SE
## H 0.04326494 0.004297766
## A 1.00000000 0.000000000# PDF page 13
for (ks in c(0.5 , 0.75 , 1)) {
result <-
svyafc(
~ wage + collgrad + hours ,
design = des_nlsw88 ,
cutoffs = list(4, 'college grad' , 26) ,
k = ks ,
g = 0 ,
dimw = c(0.5 , 0.25 , 0.25) ,
na.rm = TRUE
)
print(result)
print(attr(result, 'extra'))
}## alkire-foster SE
## [1,] 0.19135 0.0069
## coef SE
## H 0.2689563 0.009366804
## A 0.7114428 0.006847428
## alkire-foster SE
## [1,] 0.14897 0.0067
## coef SE
## H 0.1842105 0.008188892
## A 0.8087167 0.005216042
## alkire-foster SE
## [1,] 0.043265 0.0043
## coef SE
## H 0.04326494 0.004297766
## A 1.00000000 0.000000000# PDF page 17
result <-
svyafc(
~ wage + collgrad + hours,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = 1 / 3,
g = 0,
na.rm = T
)
result## alkire-foster SE
## [1,] 0.36991 0.0053## coef SE
## H 0.8082070 0.008316807
## A 0.4576895 0.004573443svyafcdec(
~ wage + collgrad + hours,
by = ~ 1,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = 1 / 3,
g = 0,
na.rm = T
)[2]## $`raw headcount ratio`
## raw headcount SE
## wage 0.19492 0.0084
## collgrad 0.76316 0.0090
## hours 0.15165 0.0076# PDF page 19
# Censored headcount ratios
for (ks in seq(1 / 3, 1, 1 / 3)) {
print (
svyafcdec(
~ wage + collgrad + hours,
by = ~ 1,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = ks,
g = 0,
na.rm = T
)[3]
)
}## $`censored headcount ratio`
## cens. headcount SE
## wage 0.19492 0.0084
## collgrad 0.76316 0.0090
## hours 0.15165 0.0076
##
## $`censored headcount ratio`
## cens. headcount SE
## wage 0.18421 0.0082
## collgrad 0.24799 0.0091
## hours 0.12756 0.0070
##
## $`censored headcount ratio`
## cens. headcount SE
## wage 0.043265 0.0043
## collgrad 0.043265 0.0043
## hours 0.043265 0.0043# Percentage contribution by dimensions
for (ks in seq(1 / 3, 1, 1 / 3)) {
print (
svyafcdec(
~ wage + collgrad + hours,
by = ~ 1,
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = ks,
g = 0,
na.rm = T
)[4]
)
}## $`percentual contribution per dimension`
## dim. % contribution SE
## wage 0.17564 0.0061
## collgrad 0.68770 0.0077
## hours 0.13666 0.0059
##
## $`percentual contribution per dimension`
## dim. % contribution SE
## wage 0.32908 0.0083
## collgrad 0.44303 0.0047
## hours 0.22789 0.0090
##
## $`percentual contribution per dimension`
## dim. % contribution SE
## wage 0.33333 0
## collgrad 0.33333 0
## hours 0.33333 0# PDF page 22
for (ks in seq(1 / 3, 1, 1 / 3)) {
print(
svyafcdec(
~ wage + collgrad + hours,
subgroup = ~ factor(married),
design = des_nlsw88,
cutoffs = list(4, 'college grad' , 26),
k = ks ,
g = 0,
dimw = NULL,
na.rm = T
)[6]
)
}## $`percentual contribution per subgroup`
## grp. % contribution SE
## 0 0.34204 0.012
## 1 0.65796 0.012
##
## $`percentual contribution per subgroup`
## grp. % contribution SE
## 0 0.32032 0.0197
## 1 0.67968 0.0197
##
## $`percentual contribution per subgroup`
## grp. % contribution SE
## 0 0.3299 0.0477
## 1 0.6701 0.0477