R/prior.R

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

#' Define a normal prior used for estimation
#' 
#' @export
#' @param mean a vector of prior means
#' @param sd a vector of prior standard deviations
#' @param defaultMean a (scalar) default mean for all elements not defined in \code{mean}
#' @param defaultSD a (scalar) default standard deviation for all elements not defined in \code{sd}
#' @details 
#' This function should never be used on its own, but only for the \code{prior} argument of 
#' \code{\link{gibbs}}. The default uses a weakly informative prior. 
#' Both arguments \code{mean} and \code{sd} are either unnamed vectors of length
#' 1, or named vectors that may contain up to three elements: \code{gamma}, \code{beta}, or 
#' \code{beta1}, and \code{beta2}. If the vector is unnamed, then the value applies to all sets of 
#' parameters. If it is named, it applies to the set of named parameters (the name \code{beta} applies 
#' to both \code{beta1} and \code{beta2} parameters). If some parameter is not defined, then the 
#' corresponding \code{default*} argument is used as a fallback. 
#' @examples
#' # default
#' normal_prior() 
#' # change the default mean and sd for all parameters
#' normal_prior(1, 2) 
#' # change the default mean for gamma, and use a different fallback mean
#' normal_prior(mean = c("gamma" = 1), defaultMean = 2) 
#' # change the default mean for both beta1 and beta2
#' normal_prior(mean = c("beta" = 2))
#' # change the default mean for all parameters
#' normal_prior(mean = c("gamma" = 1, "beta1" = 2, "beta2" = 3))
normal_prior <- function(mean = 0, sd = 5, defaultMean = 0, defaultSD = 5) {
  mean <- formatParam(mean, "mean", defaultMean)
  sd <- formatParam(sd, "sd", defaultSD)
  out <- list(mean = mean, sd = sd)
  class(out) <- "normal_prior"
  out
}

formatParam <- function(v, param, default) {
  if(length(v) > 3) stop(sprintf("Supply at most 3 %s", param))
  nams <- names(v)
  if(length(v) == 1 & is.null(nams)) {
    return(list(beta1 = v, beta2 = v, gamma = v))
  }
  if(length(nams) != length(unique(nams))) {
    stop(sprintf("Names of parameter %s are not unique", param))
  }
  if(any(! nams %in% c("gamma", "beta1", "beta2", "beta"))) {
    stop(sprintf("Names of parameter %s should be in gamma, beta1, beta2, beta", param))
  }
  if("beta" %in% nams & ("beta1" %in% nams | "beta2" %in% nams)) {
    stop(sprintf("Names of parameter %s have both beta and beta1 or beta2", param))
  }
  if("beta" %in% nams) {
    vBeta <- v["beta"]
    v <- v[nams != "beta"]
    vBeta <- rep(vBeta, 2)
    names(vBeta) <- c("beta1", "beta2")
    v <- c(v, vBeta)
    nams <- names(v)
  }
  v <- as.list(v[sort(nams)])
  for(nam in c("beta1", "beta2", "gamma")) {
    if(is.null(v[[nam]])) v[[nam]] <- default
  }
  v
}



makePrior <- function(protoPrior, kGamma, kBeta) {
  list(
    gamma = rep(protoPrior$mean$gamma, kGamma), 
    beta = cbind(
      rep(protoPrior$mean$beta1, kBeta), 
      rep(protoPrior$mean$beta2, kBeta)
    ), 
    VGamma = diag(protoPrior$sd$gamma^2, kGamma), 
    VBeta0 = diag(protoPrior$sd$beta1^2, kBeta), 
    VBeta1 = diag(protoPrior$sd$beta2^2, kBeta)
  )
}

#' @export
informative_prior <- function(formula, data, formulaMun, dataMun, munCol = j, yearCol = year) {
  modLogit <- glm(formulaMun, data = dataMun, family = binomial)
  modMultinom <- cmultinom(formula, data, munCol = munCol, yearCol = yearCol)
  vcMultinom <- vcov(modMultinom)
  k1 <- nrow(vcMultinom)/2
  k2 <- nrow(vcMultinom)
  out <- list(
    gamma = coef(modLogit), 
    beta = matrix(coef(modMultinom), ncol = 2), 
    VGamma = sandwich::vcovHC(modLogit, type = "HC1"), 
    VBeta0 = vcMultinom[1:k1,1:k1], 
    VBeta1 = vcMultinom[(k1+1):k2,(k1+1):k2]
  )
  class(out) <- "informative_prior"
  out
}