## 6.2 Alkire-Foster Class Decomposition Function Definition

svyafcdec <-
function(formula,
subgroup = ~ 1 ,
design,
g ,
cutoffs ,
k ,
dimw = NULL,
na.rm = FALSE,
...) {
if (k <= 0 |
k > 1)
stop("This functions is only defined for k in (0,1].")
if (g < 0)
stop("This function is undefined for g < 0.")
if (!is.list(cutoffs))
stop("The parameter 'cutoffs' has to be a list.")
if (length(cutoffs) != length(all.vars(formula)))
stop("number of variables in formula must exactly match cutoffs")
if (!is.null(dimw)) {
if (any(is.na(dimw))) {
stop("Invalid value in dimension weights vector.")
}
if (sum(dimw) > 1) {
stop("The sum of dimension weigths have to be equal to one.")
}
if (any(dimw > 1 |
dimw < 0)) {
stop("Dim. weights have to be within interval [0,1].")
}
}

ach.matrix <-
model.frame(formula, design$variables, na.action = na.pass)[, ] if (!is.null(dimw)) { if (any(is.na(dimw))) { stop("Invalid value in dimension weights vector.") } if (sum(dimw) > 1) { stop("The sum of dimension weigths have to be equal to one.") } if (any(dimw > 1 | dimw < 0)) { stop("Dim. weights have to be within interval [0,1].") } if (length(dimw) != ncol(ach.matrix)) { stop("Dimension weights' length differs from number of dimensions in formula") } } grpvar <- model.frame(subgroup, design$variables, na.action = na.pass)[, ]

var.class <- lapply(ach.matrix, function(x)
class(x)[1])
var.class <-
matrix(
var.class,
nrow = 1,
ncol = ncol(ach.matrix),
dimnames = list(c("var.class"), colnames(ach.matrix))
)

if (any(!(var.class %in% c("numeric", "ordered")))) {
stop(
"This function is only applicable to variables of types 'numeric' or 'ordered factor'."
)
}

w <- 1 / design$prob if (any(ach.matrix[w != 0, var.class == "numeric"] < 0, na.rm = TRUE)) stop( "The Alkire-Foster multidimensional poverty decompostition is defined for non-negative numeric variables only." ) ach.matrix <- model.frame(formula, design$variables, na.action = na.pass)[, ]
grpvar <-
model.frame(subgroup, design$variables, na.action = na.pass)[, ] if (class(grpvar) == "labelled") { stop("This function does not support 'labelled' variables. Try factor().") } if (na.rm) { nas <- apply(cbind(ach.matrix, grpvar), 1, function(x) any(is.na(x))) design <- design[nas == 0,] } w <- 1 / design$prob
ach.matrix <-
model.frame(formula, design$variables, na.action = na.pass)[, ] grpvar <- model.frame(subgroup, design$variables, na.action = na.pass)[, ]

# Deprivation Matrix
dep.matrix <-
sapply(
seq_along(var.class),
FUN = function(i) {
1 * (cutoffs[[i]] > ach.matrix[, i])
}
)
colnames(dep.matrix) <- colnames(var.class)

# Unweighted count of deprivations:
# depr.count <- rowSums( dep.matrix )

# deprivation k cut
if (is.null(dimw)) {
dimw = rep(1 / ncol(dep.matrix), ncol(dep.matrix))
}

# Weighted sum of deprivations:
depr.sums <- rowSums(sweep(dep.matrix, MARGIN = 2 , dimw, *))

# k multidimensional cutoff:
multi.cut <- depr.sums * (depr.sums >= k)
#rm(dep.matrix)

# Censored Deprivation Matrix
cen.dep.matrix <-
sapply(
seq_along(cutoffs) ,
FUN = function(x) {
if (var.class[x] == "numeric") {
1 * (cutoffs[[x]] > ach.matrix[, x]) * ((cutoffs[[x]] - ach.matrix[, x]) / cutoffs[[x]]) ^
g
} else {
1 * (cutoffs[[x]] > ach.matrix[, x])
}
}
)
colnames(cen.dep.matrix) <- colnames(var.class)

cen.dep.matrix[multi.cut == 0 & !is.na(multi.cut),] <- 0

# Sum of censored deprivations:
cen.depr.sums <-
rowSums(sweep(cen.dep.matrix, MARGIN = 2 , dimw, *))

if (any(is.na(cen.depr.sums)[w > 0])) {
# overall result:
overall.result <-
matrix(
c(NA, NA),
nrow = 1,
ncol = 2,
dimnames = list("overall", c("alkire-foster", "SE"))
)

# group breakdown:
if (!is.null(levels(grpvar))) {
grp.estim <-
matrix(
rep(NA, 2 * length(levels(grpvar))),
nrow = length(levels(grpvar)),
ncol = 2,
dimnames = list(levels(grpvar), c("alkire-foster", "SE"))
)
grp.contr <-
matrix(
rep(NA, 2 * length(levels(grpvar))),
nrow = length(levels(grpvar)),
ncol = 2,
dimnames = list(levels(grpvar), c("contribution", "SE"))
)

}

# dimensional breakdown:
dim.raw.hc <-
matrix(
rep(NA, 2 * ncol(cen.dep.matrix)),
nrow = ncol(cen.dep.matrix),
ncol = 2,
dimnames = list(colnames(cen.dep.matrix), c("raw headcount", "SE"))
)
dim.cen.hc <-
matrix(
rep(NA, 2 * ncol(cen.dep.matrix)),
nrow = ncol(cen.dep.matrix),
ncol = 2,
dimnames = list(colnames(cen.dep.matrix), c("cens. headcount", "SE"))
)
dim.contri <-
matrix(
rep(NA, 2 * ncol(cen.dep.matrix)),
nrow = ncol(cen.dep.matrix),
ncol = 2,
dimnames = list(colnames(cen.dep.matrix), c("contribution", "SE"))
)

# set up result object:
if (is.null(levels(grpvar))) {
rval <- list(overall.result, dim.raw.hc, dim.cen.hc, dim.contri)
} else {
rval <-
list(overall.result,
dim.raw.hc,
dim.cen.hc,
dim.contri,
grp.estim,
grp.contr)
}

return(rval)

}

raw.hc.ratio <-
survey::svymean(dep.matrix[, ] , design, na.rm = TRUE)
cen.hc.ratio <-
survey::svymean(cen.dep.matrix[, ] , design, na.rm = TRUE)

U_0 <-
list(value = sum(w[w > 0]), lin = rep(1, length(w)))
U_1 <-
list(value = sum(w[w > 0] * cen.depr.sums[w > 0]), lin = cen.depr.sums)

# overall alkire-foster index:
overall <- survey::svymean(cen.depr.sums, design, na.rm = TRUE)
names(overall)[1] <-
attr(overall, "statistic") <- "alkire-foster"

# group decomposition
if (!is.null(levels(grpvar))) {
grp.pctg.estm <- NULL
grp.pctg.estm_var <- NULL
grp.pctg.cont <- NULL
grp.pctg.cont_lin <-
matrix(data = rep(w, length(levels(grpvar))), nrow = length(w))
for (i in seq_along(levels(grpvar))) {
w_i <- w * (grpvar == levels(grpvar)[i])
U_1_i <-
list(value = sum(cen.depr.sums[w_i > 0] * w_i[w_i > 0]),
lin = cen.depr.sums * (grpvar == levels(grpvar)[i]))
U_0_i <-
list(value = sum(w_i[w_i > 0]),
lin = 1 * (grpvar == levels(grpvar)[i]))

list_all <-
list(
U_0 = U_0,
U_1 = U_1,
U_0_i = U_0_i,
U_1_i = U_1_i
)

estimate <-
contrastinf(quote((U_1_i / U_0_i)), list_all)
grp.pctg.estm[i] <- estimate$value estimate$lin[is.na(estimate$lin)] <- 0 grp.pctg.estm_var[i] <- survey::svyrecvar( estimate$lin * w_i,
design$cluster, design$strata,
design$fpc, postStrata = design$postStrata
)

estimate <-
contrastinf(quote((U_0_i / U_0) * (U_1_i / U_0_i) / (U_1 / U_0)), list_all)
grp.pctg.cont[i] <- estimate$value estimate$lin[is.na(estimate$lin)] <- 0 grp.pctg.cont_lin[, i] <- estimate$lin

}

grp.pctg.cont_var <-
survey::svyrecvar(
grp.pctg.cont_lin * w,
design$cluster, design$strata,
design$fpc, postStrata = design$postStrata
)

grp.contr.estimate <-
matrix(grp.pctg.estm,
ncol = 1,
dimnames = list(levels(grpvar), "alkire-foster"))
attr(grp.contr.estimate, "names") <- levels(grpvar)
attr(grp.contr.estimate, "var") <- grp.pctg.estm_var
attr(grp.contr.estimate, "statistic") <- "alkire-foster"
class(grp.contr.estimate) <- c("svystat")

grp.contr.pct <-
matrix(
grp.pctg.cont,
ncol = 1,
dimnames = list(levels(grpvar), "grp. % contribution")
)
attr(grp.contr.pct, "names") <- levels(grpvar)
attr(grp.contr.pct, "var") <- grp.pctg.cont_var
attr(grp.contr.pct, "statistic") <- "grp. % contribution"
class(grp.contr.pct) <- c("svystat")

rm(
grp.pctg.estm,
grp.pctg.estm_var,
grp.pctg.cont,
grp.pctg.cont_lin,
grp.pctg.cont_var
)

}

# dimensional decomposition:
dim.contr <- NULL
dim.contr_lin <-
matrix(data = rep(w, ncol(cen.dep.matrix)), nrow = length(w))
for (i in 1:ncol(cen.dep.matrix)) {
wj <-
list(value = dimw[i], lin = rep(0, length(cen.depr.sums)))

H_0 <-
list(value = sum(cen.dep.matrix[w > 0 , i] * w[w > 0]),
lin = cen.dep.matrix[, i] * (w > 0))
list_all <- list(
U_0 = U_0,
U_1 = U_1,
H_0 = H_0,
wj = wj
)
estimate <-
contrastinf(quote(wj * (H_0 / U_0) / (U_1 / U_0)), list_all)

dim.contr[i] <- estimate$value estimate$lin[is.na(estimate$lin)] <- 0 dim.contr_lin[, i] <- estimate$lin

}

dim.contr_var <-
survey::svyrecvar(
dim.contr_lin * w,
design$cluster, design$strata,
design$fpc, postStrata = design$postStrata
)

dim.result <- dim.contr
attr(dim.result, "names") <- colnames(ach.matrix)
attr(dim.result, "var") <- dim.contr_var
attr(dim.result, "statistic") <- "dim. % contribution"
class(dim.result) <- c("svystat")

rm(dim.contr, dim.contr_lin, dim.contr_var)

# set up result object:
if (is.null(levels(grpvar))) {
rval <- list(overall, raw.hc.ratio, cen.hc.ratio, dim.result)
names(rval) <-
list(
"overall",
"percentual contribution per dimension"
)

} else {
rval <-
list(
overall,
raw.hc.ratio,
cen.hc.ratio,
dim.result,
grp.contr.estimate,
grp.contr.pct
)
names(rval) <-
list(
"overall",
}