R/fm_ZBM.R

In FuzzyPovertyR: Estimation of Fuzzy Poverty Measures

# Fuzzy monetary poverty estimation
#
# @description Calculates bootstrap percentiles from Zedini and Belhadj (2015)
#
# @param x A numeric vector of a predicate variable (or poverty predicate)
# @param R The number of bootstrap replicates (defaults to 500)
#
# @return A matrix of bootstrapped percentiles
#
# @references
# Zedini, A., & Belhadj, B. (2015). A New Approach to Unidimensional Poverty Analysis: Application to the T unisian Case. Review of Income and Wealth, 61(3), 465-476.

bootP <- function(x, R = 500){
  M <- length(x)
  bootReps <- replicate(R, sample(x, size = M, replace = TRUE))
  boot.p <- apply(bootReps, 2, quantile, probs = seq(0.01, 1, 0.01))
  P <- apply(boot.p, 1, mean)
  return(P)
}

# Fuzzy monetary poverty estimation
#
# @param x A numeric vector of a predicate variable (or poverty predicate)
# @param a Membership function parameter
# @param b Membership function parameter
# @param c Membership function parameter
#
# @return The membership grade
#
# @references
# Zedini, A., & Belhadj, B. (2015). A New Approach to Unidimensional Poverty Analysis: Application to the Tunisian Case. Review of Income and Wealth, 61(3), 465-476.

FN <- function(x, a, b, c){
  y <- rep(NA_real_, length(x))
  y[x<a | x >= c] <- 0
  y[a<=x & x<= b] <- 1
  y[b<=x & x < c] = x[b<=x & x<c]*(-1/(c-b))+(c/(c-b))
  return(y)
}

# Fuzzy monetary poverty estimation
#
# @param x A numeric vector of a predicate variable (or poverty predicate)
# @param P A matrix of bootstapped percentiles
#
# @return The membership grade matrix

MemberhsipGradesMatrix <- function(x, P){
  n <- length(x)
  MGM <- matrix(0, nrow = n, ncol = 1)
  MGM[,1] <- FN(x = x, a = 0, b = P[1], c = P[2])
  MGM <- cbind(MGM, do.call(cbind, lapply(2:99, function(j) FN(x = x, a = P[j-1], b = P[j], c = P[j+1]))))
  MGM <- cbind(MGM, FN(x = x, a = P[99], b = 0.5*(P[99] + P[100]), c = P[100]))
  return(MGM)
}

# Fuzzy monetary poverty estimation
#
# @param x A numeric vector of fuzzy numbers
#
# @return the convolution of fuzzy numbers

Fuzzy_conv <- function(x) {
  Min <- Vectorize(function(a,b) min(a,b), vectorize.args = "b")
  Max <- Vectorize(function(a,b) max(a,b), vectorize.args = "b")
  l <- length(x)
  m <- Min(x[1], x[2:l])
  P <- Max(m[1], m[2:length(m)])
  return(P)
}

# Fuzzy monetary poverty estimation
#
# @param predicate A numeric vector of a predicate variable (or poverty predicate)
# @param hh.size A numeric vector with the size of the household
# @param weight A numeric vector of sampling weights. if NULL simple random sampling weights will be used
# @param breakdown A factor of sub-domains to calculate estimates for (using the same alpha).
# @param ID A numeric or character vector of IDs. if NULL (the default) it is set as the row sequence.
#
# @import ecp
# @return The membership grades for each poverty state
#
# @references
# Zedini, A., & Belhadj, B. (2015). A New Approach to Unidimensional Poverty Analysis: Application to the Tunisian Case. Review of Income and Wealth, 61(3), 465-476.
# Belhadj, B., & Matoussi, M. S. (2010). Poverty in tunisia: A fuzzy measurement approach. Swiss Journal of Economics and Statistics, 146(2), 431-450.

fm_ZBM <- function(predicate, hh.size, weight, breakdown , ID){
  # Zendini, Belhadi, Matoussi
  k <- 3 # The number of change points locations to estimate
  N <- length(predicate)
  if (is.null(ID)) ID <- seq_len(N)
  if(is.null(hh.size)) hh.size <- rep(1, N)

  #--- Step 0 ---#
  P <- bootP(x = predicate) # bootstrap estimates of percentiles
  a = P[99]; b = 0.5*(P[99] + P[100]); c = P[100]
  MGM <- MemberhsipGradesMatrix(x = predicate, P = P)

  #--- Step 1 ---#
  cluster.idx <- lapply(1:N, function(i) ecp::e.divisive( X = matrix(MGM[i,], ncol = 1), k = k)$cluster)

  #--- Step 2-3 ---#
  step2_3 <- lapply(1:N, function(i) tapply( MGM[i,], cluster.idx[[i]], max))
  step2_3.mat <- matrix(NA_real_, nrow = N, ncol = k, byrow = T, dimnames = list(NULL, 1:k))
  for(i in 1:N) step2_3.mat[i, names(step2_3[[i]])] = step2_3[[i]]
  # if(sum(!complete.cases(step2_3.mat))!=0) warning(paste("for", sum(!complete.cases(step2_3.mat)), "cases NA were generated by e.divisive\n"))
  step2_3.mat[is.na(step2_3.mat)] = 0

  #--- Step 4 ---#

  if(!is.null(breakdown)) {

    states.split <- split(data.frame(step2_3.mat, hh.size, weight, breakdown), f= ~ breakdown)
    Q <- t(sapply(states.split, function(Z) apply(Z[,1:k], 2, function(z) weighted.mean(x = z*Z[["hh.size"]], w = Z[["weight"]]))))
    colnames(Q) = c(1:k)
    P <- apply(Q, 1, Fuzzy_conv)
  } else {
        Q <- apply(step2_3.mat, 2, function(y) weighted.mean(x = y*hh.size, w = weight))
        P <- Fuzzy_conv(Q)
  }
  fm_data <- data.frame(ID = ID, predicate = predicate, weight = weight)
  order.idx <- order(fm_data$predicate)
  fm_data$mu = FN(fm_data$predicate,
                  a=a,
                  b=b,
                  c=c)
  fm_data <- fm_data[order.idx,]
  return(list( results = fm_data,
               MGM = MGM,
               mu_k =step2_3.mat[order.idx,],
               estimate = Q, parameters = list(P = P, k = k, a  = a, b = b, c = c), fm = "ZBM"))
}