R/helpers.R

In rferrali/rogali: A series of functions for the corrupt_audit project

getNObs <- function(formula, data, munCol, yearCol) {
  munCol <- rlang::enquo(munCol)
  yearCol <- rlang::enquo(yearCol)
  y <- model.response(model.frame(formula, data))
  nSpell <- sum(y)
  
  tmp <- data %>% 
    dplyr::rename(mun = !!munCol, year = !!yearCol) %>% 
    dplyr::mutate(year = as.numeric(as.character(year)))
  
  nMun <- length(unique(tmp$mun))
  
  nObs <- tmp %>% 
    dplyr::group_by(mun) %>% 
    dplyr::mutate(keep = year == min(year)) %>% 
    dplyr::pull(keep)
  nObs <- sum(y[nObs,])
  
  list(nSpell = nSpell, nObs = nObs, nMun = nMun)
}

poolDf <- function(formula, df, weights = NULL) {
  y <- model.response(model.frame(formula, df)) 
  if(!is.null(weights)) {
    y <- y * weights
    keep <- weights > 0
    y <- y[keep,,drop = F]
    df <- df[keep,]
  }
  y <- as.data.frame(y)
  df <- model.frame(delete.response(terms(formula)), df)
  df <- cbind(y, df)
  df %>% 
    dplyr::group_by_at(-(1:3)) %>% 
    dplyr::summarize_all(sum) %>% 
    dplyr::ungroup()
}


getBoots <- function(bootSrc, seeds) {
  sapply(seeds, function(seed, bootSrc) {
    set.seed(seed)
    unlist(tapply(bootSrc$mun, bootSrc$state, function(v) sample(v, length(v), replace = T)))
  }, bootSrc = bootSrc, simplify = F)
}

linCoefs <- function(mod) {
  coefs <- as.numeric(t(coef(mod)))
  coefNames <- colnames(coef(mod))
  coefNames <- c(
    paste0("beta1-", coefNames), 
    paste0("beta2-", coefNames)
  )
  names(coefs) <- coefNames
  coefs
}

bootCoefAdjust <- function(coefs, coefNames) {
  sto <- rep(NA, length(coefNames))
  m <- match(names(coefs), coefNames)
  sto[m] <- coefs
  names(sto) <- coefNames
  sto
}



se <- function(object, ... = ...) {
  if(is.null(vcov(object))) {
    return(NULL)
  } else {
    se <- sqrt(diag(vcov(object)))
    names(se) <- names(coef(object))
    if(...length() > 0) {
      se <- tibble::as_tibble(t(se))
      se <- dplyr::select(se, ... = ...)
      cnames <- colnames(se)
      se <- as.numeric(se)
      names(se) <- cnames
    }
    se
  }
}

#' Posterior probabilities Pr(Z = 1)
#' @param object a \code{gibbs} object
#' @param burn number or percent of iterations to discard
#' @param thin thinning interval. Keep every \code{thin}-th iteration. 
#' @param cutoff cutoff in probability for Z = 1
#' @param ... additional arguments to be passed to or from methods
#' @return \code{posterior_prob} returns a \code{tibble} with the observed 
#' type and posterior probability of being of type 1. \code{confusion_table} 
#' returns a \code{table} comparing predicted probabilities and types for 
#' those municipalities whose type is observed using the cutoff set in 
#' \code{cutoff}. 
#' @name posteriorprobs
NULL

#' @rdname posteriorprobs
#' @export
posterior_probs <- function(object, ...) {
  UseMethod("posterior_probs", object)
}

#' Confusion table for posterior probabilities
#' @rdname posteriorprobs
#' @export
confusion_table <- function(object, cutoff = .5, ...) {
  UseMethod("confusion_table", object)
}