R/mid.R

In MLID: Multilevel Index of Dissimilarity

## These are the functions for fitting the multilevel index using lme4,
## calculating the percentage of the variance at each level and
## also the holdback scores


.mid <- function(data, vars = c(1, 2), levels, expected = FALSE, nsims = 100) {
  if (is.character((levels))) {
    ifelse (all(levels %in% names(data)), levels <- match(levels, names(data)),
            stop("Higher level grouping variable not found"))
  }
  if (anyNA(data[,levels])) stop("Data contain NAs")
  if (any(sapply(data[, levels], is.numeric)))
    warning("The levels contain numeric data")
  id <- .idx(data, vars, expected, nsims)
  ols <- attr(id, "ols")
  vv <- attr(id, "vars")
  df <- attr(id, "data")
  lvls <- names(data)[levels]
  data <- data.frame(y = data[, vars[1]], x = data[, vars[2]], data[, levels])
  data$y <- data$y / sum(data$y)
  data$x <- data$x / sum(data$x)
  f <- "y ~ 0"
  for(k in 3: ncol(data)) {
    f <- paste(f, "+", paste0("(1|", names(data)[k],")"))
  }
  mlm <- lme4::lmer(f, data=data, offset=data$x)
  attributes(id) <- list(ols = ols, vars = vv, levels = lvls,
                         mlm = mlm, variance = .varshare(mlm),
                         holdback = holdback(mlm),
                         data = df)
  class(id) <- "index"
  return(id)
}


.mlvar <- function(mlm) {

  vv <- lme4::VarCorr(mlm)
  res <- attr(vv, "sc")

  for(i in 1: length(vv)) {

    res <- c(res,attr(vv[[i]], "stddev"))

  }

  res <- res^2
  names(res) <- c("Base",names(vv))
  return(res)
}


.varshare <- function(mlm) {

  mlvar <- .mlvar(mlm)
  mlvar <- mlvar/sum(mlvar)*100
  return(round(mlvar,2))

}


#' Holdback scores
#'
#' Calculates the holdback scores for a multilevel index of dissimilarity
#'
#' For the index of dissimilarity (ID), the residuals are the differences
#' between the share of the Y population and the share of the X population per
#' neighbourhood. For the multilevel index, the residuals are estimated at and
#' partitioned between each level of the model. The holdback scores take each
#' level in the model in turn and set the residuals (the effects) at that
#' level to zero, then recalculating the ID on that basis and recording
#' the percentage change in the original value that occurs. The holdback
#' scores are calculated automatically as part of the function \code{\link{id}}
#' and can be viewed through \code{print(index)}, where \code{index} is the
#' object returned by the function, or as \code{attr(index, "holdback")}.
#'
#' @param mlm an object of class \code{lmerMod} generated by the
#' \code{lme4} package
#' @return a numeric vector containing the holdback scores
#' @seealso \code{\link{id}} \code{\link{print.index}}


holdback <- function(mlm) {
  rvals <- .rvals(mlm)
  k <- ncol(rvals)
  hb <- rep(NA, k)
  rwsm <- rowSums(rvals)
  ID <- sum(abs(rwsm))
  for(i in 1:k) {

    rwsm <- rowSums(rvals[,-i])
    hb[i] <-  sum(abs(rwsm))/ID - 1

  }
  names(hb) <- colnames(rvals)
  return(round(hb*100,2))

}