R/residuals.bayesbr.R

In bayesbr: Beta Regression on a Bayesian Model

#' @title Residuals for \code{bayesbr} Objects
#' @aliases residuals.bayesbr
#' @name residuals.bayesbr
#' @description A function that receives model information and calculates the residuals according to the required residual.
#' @usage \method{residuals}{bayesbr}(object, type = c("", "quantile", "sweighted", "pearson","ordinary"),...)
#' @param object an object of the class \emph{bayesbr}, containing the list returned from the \code{\link{bayesbr}} function.
#' @param type A character containing the residual type returned by the model among the possibilities. The type of residue can be \emph{quantile}, \emph{sweighted}, \emph{pearson} or \emph{ordinary}. The default is \emph{quantile}.
#'@param ... further arguments passed to or from other methods.
#' @details
#'The definitions of the waste generated by the package are available in Espinheira (2008): "pearson" in Equation 2, "sweighted" in Equation 7; and in Pereira (2019): "quantile" in Equation 5;
#'
#'The type of residue "response" is calculated from the difference between the estimated theta and the variable response of the model.
#'@return A vector containing the model residual according to the type of residual calculated
#'@references
#'\doi{10.1080/0266476042000214501} Ferrari, S., & Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. \emph{Journal of applied statistics}, \bold{31}(7), 799-815.
#'
#'\doi{10.1080/00949650701829380} Simas, A. B., & Cordeiro, G. M. (2009). Adjusted Pearson residuals in exponential family nonlinear models. \emph{Journal of Statistical Computation and Simulation}, \bold{79}(4), 411-425.
#'
#'\doi{10.1080/02664760701834931} Espinheira, P. L., Ferrari, S. L., & Cribari-Neto, F. (2008). On beta regression residuals. \emph{Journal of Applied Statistics}, \bold{35}(4), 407-419.
#'
#'\doi{10.1080/00949655.2012.736993} Anholeto, T., Sandoval, M. C., & Botter, D. A. (2014). Adjusted Pearson residuals in beta regression models. \emph{Journal of Statistical Computation and Simulation}, \bold{84}(5), 999-1014.
#'
#'\doi{10.1080/03610918.2017.1381740} Pereira, G. H. (2019). On quantile residuals in beta regression. \emph{Communications in Statistics-Simulation and Computation}, \bold{48}(1), 302-316.
#'@seealso \code{\link{bayesbr}},\code{\link{summary.bayesbr}},\code{\link{predict.bayesbr}}
#'@examples
#'data("CarTask", package = "bayesbr")
#'
#'bbr = bayesbr(probability~task + NFCCscale,data=CarTask,
#'              iter = 100, mean_betas = c(1, 0.5,1.2))
#'
#'residuals(bbr, type = "quantile")
#'residuals(bbr, type = "ordinary")
#'residuals(bbr, type = "sweighted")
#'residuals(bbr, type = "pearson")
#'@export
residuals.bayesbr = function(object,type = c("","quantile","sweighted", "pearson","ordinary"),...){
  type = match.arg(type)
  if(type==""){
    type = object$residuals.type
  }
  Y = object$y
  p = object$info$p
  q = object$info$q
  n = object$info$n
  theta = object$info$samples$theta
  zeta = object$info$samples$zeta
  res = NULL
  if(type=="quantile"){
    sum = 0
    res_vector = c()
    for (i in 1:n) {
      v_theta = c()
      v_zeta = c()
      if(p==0){
        v_theta = theta$theta
      }
      else{
        aux = paste0('theta[',i,']')
        v_theta = theta[[aux]]
      }
      if(q==0){
        v_zeta = zeta$zeta
      }
      else{
        aux = paste0('zeta[',i,']')
        v_zeta = zeta[[aux]]
      }
      mean_theta = mean(v_theta)
      mean_zeta = mean(v_zeta)
      A = mean_theta*mean_zeta
      B = mean_zeta-A
      res = qnorm(pbeta(q = Y[i], shape1 = A,shape2 = B))
      res_vector = c(res_vector,res)
    }
  }
  if(type=="sweighted"){
    res_vector = c()
    for (i in 1:n) {
      v_theta = c()
      v_zeta = c()
      if(p==0){
        v_theta = theta$theta
      }
      else{
        aux = paste0('theta[',i,']')
        v_theta = theta[[aux]]
      }
      if(q==0){
        v_zeta = zeta$zeta
      }
      else{
        aux = paste0('zeta[',i,']')
        v_zeta = zeta[[aux]]
      }
      mean_theta = mean(v_theta)
      mean_zeta = mean(v_zeta)
      y_star = log(Y[i]/(1-Y[i]))
      mean_star = (digamma(mean_theta*mean_zeta) - digamma((1-mean_theta)*mean_zeta) )
      variance_star = (trigamma(mean_theta*mean_zeta) + trigamma((1-mean_theta)*mean_zeta) )
      res = (y_star-mean_star)/(variance_star^(1/2))
      res_vector = c(res_vector,res)
    }
  }

  #if(type=="sweighted2") a fazer

  if(type=="pearson"){
    res_vector = c()
    for (i in 1:n) {
      v_theta = c()
      v_zeta = c()
      if(p==0){
        v_theta = theta$theta
      }
      else{
        aux = paste0('theta[',i,']')
        v_theta = theta[[aux]]
      }
      if(q==0){
        v_zeta = zeta$zeta
      }
      else{
        aux = paste0('zeta[',i,']')
        v_zeta = zeta[[aux]]
      }
      mean_theta = mean(v_theta)
      mean_zeta = mean(v_zeta)
      diff = (Y[i] - mean_theta)
      variance = (mean_theta*(1-mean_theta))/(1+mean_zeta)
      res = diff/(variance^(1/2))
      res_vector = c(res_vector,res)
    }
  }
  if(type=="ordinary"){
    res_vector = c()
    for (i in 1:n) {
      v_theta = c()
      v_zeta = c()
      if(p==0){
        v_theta = theta$theta
      }
      else{
        aux = paste0('theta[',i,']')
        v_theta = theta[[aux]]
      }
      if(q==0){
        v_zeta = zeta$zeta
      }
      else{
        aux = paste0('zeta[',i,']')
        v_zeta = zeta[[aux]]
      }
      mean_theta = mean(v_theta)
      mean_zeta = mean(v_zeta)
      diff = (Y[i] - mean_theta)
      res_vector = c(res_vector,diff)
    }
  }
  names(res_vector) = 1:n
  return(res_vector)
}