R/helpers.R

In lfpdroubi/appraiseR: R Package For Real Estate Appraisals

#' @import stats
#' @import ggplot2

precision = getFromNamespace("precision", "scales")

#' @importFrom magrittr %>%
#' @export
magrittr::`%>%`

#' @importFrom magrittr %<>%
#' @export
magrittr::`%<>%`

#' Parent Function: Power of a number.
#'
#' This is an internal function to generate another ones which will be
#' effectivelly exported. This function is not exported in NAMESPACE, therefore
#' it is not available for the end-user.
#'
#' @param exponent The exponent.
#' @return A parent function that allows the user to create a closure that
#'   returns the n-th power of its argument.
#' @examples
#'
#' \dontrun{
#' power <- function(exponent) {
#' function(x) {
#'   x ^ exponent
#' }
#' }
#'
#' square <- power(2)
#' square_root <- power(.5)
#'
#' square(2) #4
#' square_root(4) #2
#' }

power <- function(exponent) {
  function(x) {
    x^exponent
  }
}

#' Reciprocal of the square of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' rsqr(2)
#' all.equal(rsqr(2), 1/2^2)
#' rsqr(1:10)

#rsqr <- power(-2)
rsqr <- function(x) 100*x^(-2)

#' Reciprocal (1/x) of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' rec(2)
#' all.equal(rec(2), 1/2)
#' rec(1:10)

#rec <- power(-1)
rec <- function(x) -100*x^(-1)

#' Reciprocal of the square root of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' rsqrt(4)
#' all.equal(rsqrt(2), 1/sqrt(2))
#' rsqrt(rsqr(1:10))

#rsqrt <- power(-.5)
rsqrt <- function(x) -10*x^(-1/2)

#' Square of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' sqr(2)
#' sqr(1:10)

sqr <- power(2)

#' Cube of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' cube(2)

cube <- power(3)

#' Cubic root of a number
#'
#' @param x a numeric vector or array
#' @export
#' @examples
#'
#' cubroot(8)

cubroot <- power(1/3)

#' Returns the inverse of the function f (character)
#'
#' @param FUN a function of the box-cox family (rsqr(), rec(), rsqrt(), log(),
#' cubroot(), sqrt(), I() and sqr())
#' @return the inverse function, in character format.
#' @export
#' @examples
#' invFunc("log")

invFunc <- function(FUN){
  switch(FUN,
         rsqr = "rsqrt",
         rec = "rec",
         rsqrt = "rsqr",
         log = "exp",
         exp = "log",
         sqrt = "sqr",
         sqr = "sqrt",
         cube = "cubroot",
         cubroot = "cube"
  )
}

#' Returns a vector generated with the inverse of the function f
#'
#' @param x A vector or object of type
#' @param FUN a function of the box-cox family (rsqr(), rec(), rsqrt(), log(),
#' cubroot(), sqrt(), I() and sqr())
#' @return a
#' @export
#' @examples
#' inverse(rsqr(10), "rsqr")
#' inverse(rec(10), "rec")
#' inverse(rsqrt(10), "rsqrt")
#' inverse(log(1), "log")
#' inverse(sqrt(4), "sqrt")
#' inverse(sqr(4), "sqr")
#' dados <- st_drop_geometry(centro_2015)
#' fit <- lm(log(valor) ~ ., data = dados)
#' aval <- new_data(fit)
#' Y <- predict(fit, newdata = aval, interval = "confidence")
#' inverse(Y, "log")
#'

inverse <- function(x, FUN) {
  switch(FUN,
         rsqr = appraiseR::rsqrt(x),
         rec = appraiseR::rec(x),
         rsqrt = appraiseR::rsqr(x),
         log = exp(x),
         exp = log(x),
         cubroot = appraiseR::cube(x),
         sqrt = appraiseR::sqr(x),
         identity = identity(x),
         sqr = sqrt(x),
         cube = appraiseR::cubroot(x)
  )
}

#' Return the median value or the modal value of a vector, depending on whether
#' the vector is numeric or a factor.
#'
#' @inheritParams stats::median
#' @return the median value for objects of the class integer or double. The
#'   modal value for objects of class factor.
#' @name centre
#' @export
centre <- function(x, ...) UseMethod("centre")

#' @rdname centre
#' @examples
#' x <- c(-3, -2, 0, NA, 1, 1, 3)
#' centre(x) # 0 (default is mean(x))
#' centre(x, FUN = median) # 0.5
#' centre(x, FUN = raster::modal) # 1
#'
#' data(centro_2015)
#' lapply(centro_2015, centre)
#' # lapply(centro_2015, FUN = centre, median)
#' lapply(centro_2015[, -7], FUN = centre, median)
#' @export
centre.numeric <- function(x, FUN = mean, na.rm = TRUE, ...) {
  FUN(x, na.rm = na.rm, ...)
}

#' @rdname centre
#' @examples
#' y <- c(-3, -2, 0, 1, 1, 3)
#' y <- as.factor(y)
#' centre(y) # the reference level
#' centre(y, FUN = raster::modal) # the most frequent (modal) value
#' library(sf)
#' data(centro_2015)
#' centre(centro_2015$padrao) # the reference level
#' centre(centro_2015$padrao, FUN = raster::modal) # the most frequent value
#' centre(centro_2015$padrao, FUN = median)
#' @export
centre.factor <- function(x, FUN, na.rm = TRUE, ...){

  if (missing(FUN)) {
    x <- levels(x)[1]
  } else {
    FUN <- match.fun(FUN)
    x <- FUN(x, na.rm = na.rm, ...)
  }
  return(x)
}
#' centre(centro_2015$geometry)
#' @export
centre.sfc_POINT <- function(x, FUN = mean, ...){
  z <- sf::st_coordinates(x)
  apply(z, 2, FUN, ...)
}

#' Extract object parameters
#'
#' Returns the parameters used to build a model.
#'
#' @param object A model object.
#' @param \dots not used.
#' @return the parameters, predictors and response names besides the
#' original data used to build the model.
#' @name parameters
#' @export
parameters <- function(object, ...) {
  UseMethod("parameters")
}

#' @rdname parameters
#' @examples
#' # 1. Real dataset
#'
#' library(sf)
#' data(centro_2015)
#' centro_2015 <- st_drop_geometry(centro_2015)
#' fit <- lm(log(valor) ~ ., centro_2015, subset = -c(31, 39))
#' p <- parameters(fit)
#' p$parameters
#' p$rhs
#' p$predictors
#' p$response
#' p$lhs
#' p$depvarTrans
#' p$data
#'
#' # 2. Null model
#'
#' fit1 <- lm(log(valor) ~ 1, data = centro_2015)
#' p1 <- parameters(fit1)
#' p1$parameters
#' p1$rhs
#' p1$predictors
#' p1$response
#' p1$lhs
#' p1$depvarTrans
#' p1$data
#'
#' # 3. Random generated data
#'
#' n <- 20
#' set.seed(1)
#' Frente <- rnorm(n = n, mean = 12, sd = 2.5)
#' Prof <- rnorm(n = n, mean = 25, sd = 5)
#' Area <- Frente*Prof
#' # quadratic relationship
#' PU <- 5000 - 10*Area + .005*Area^2 + 10*Frente + rnorm(n, mean = 0, sd = 150)
#' d <- data.frame(PU, Area, Frente, Prof)
#' fit2 <- lm(PU ~ poly(Area, 2) + Frente, data = d)
#' p2 <- parameters(fit2)
#' p2$parameters
#' p2$rhs
#' p2$predictors
#' p2$response
#' p2$lhs
#' p2$depvarTrans
#' p2$data
#'
#' fit3 <- lm(PU ~ poly(Area, 2, raw=TRUE) + Frente, data = d)
#' p3 <- parameters(fit3)
#' p3$parameters
#' p3$rhs
#' p3$predictors
#' p3$response
#' p3$lhs
#' p3$depvarTrans
#' p3$data
#' @export
#'
parameters.lm <- function(object, ...) {
  z <- object

  cl <- stats::getCall(z)
  myformula <- stats::formula(z)
  #data <- as.data.frame(eval(z$call$data))
  data <- eval(stats::getCall(z)$data, environment(formula(z)))
  vars <- all.vars(myformula)
  tt <- terms(myformula)

  resp <- all.vars(update(myformula, . ~ 1))
  preds <- setdiff(vars, resp)

  #lhs <- rownames(attr(z$terms, "factors"))[1]
  lhs <- as.character(myformula)[2]
  rhs <- as.character(myformula)[3]

  if (stringr::str_detect(lhs, paste("\\(", resp, "\\)", sep = ""))) {
    depvarTrans <- stringr::str_replace(lhs,
                                        pattern = paste("\\(", resp, "\\)",
                                                        sep = ""),
                                        replacement = "")
  } else {
    depvarTrans <- NA
  }


  param <-
    list(parameters = c(resp, preds),
         predictors = preds,
         response = resp,
         depvarTrans = depvarTrans,
         lhs = lhs,
         rhs = rhs,
         data = data,
         call = cl)

  return(param)
}

#' @rdname parameters
#' @examples
#' dados <- st_drop_geometry(centro_2015)
#' best_fit <- bestfit(valor ~ ., dados)
#' parameters(best_fit)
#' @export
#'
parameters.bestfit <- function(object, ...) {
  z <- object
  cl <- z$call
  data <- eval(cl$data, environment(stats::formula(z)))

  resp <- z$response
  preds <- z$predictors

  param <-
    list(parameters = c(resp, preds),
         predictors = preds,
         response = resp,
         data = data,
         call = cl)

  return(param)
}

#' Builds \code{newdata} argument to be used in \link{predict.lm}
#'
#' Builds a new \code{data.frame} containing only elements
#' to be appraised from the current \code{data.frame}
#'
#' @param object object of class \code{lm}
#'
#' @examples
#' dados <- st_drop_geometry(centro_2015)
#' fit <- lm(log(valor) ~ ., data = dados)
#' new_data(fit)
#' @export

new_data <- function(object) {
  z <- object
  params <- parameters(z)
  response <- params$response
  response <- as.name(response)
  parameters <- params$parameters
  data <- params$data
  aval <-
    data %>%
    dplyr::filter(is.na(!!response)) %>%
    dplyr::select(parameters)
  aval
}

#' Wrapper around format
#'
#' This function is only a wrapper around \link{format} function that uses standard
#' brazillian formats by default
#'
#' @param x a number to be formatted by \code{format}
#' @inheritParams base::format
#' @export
brformat <- function(x, decimal.mark = ",", big.mark = ".", digits = 2,
                     nsmall = 2, scientific = FALSE, ...) {
  format(x, decimal.mark = decimal.mark, big.mark = big.mark, digits = digits,
         nsmall = nsmall, scientific = scientific, ...)
}

#' @rdname brformat
#' @export
brf <- brformat

#' Wrapper around brformat
#'
#' This is a wrapper around \link{brformat}.
#'
#' @param x a real number
#' @param prefix currency units. Defaults for brazilian reais.
#' @param \ldots further arguments to be passed to \link{brformat}.
#' @return text in Reais.
#' @examples
#' Reais(100)
#' @export
#'
Reais <- function(x, prefix = "R$ ", ...) {
  paste(prefix, brformat(x, ...), sep = "")
}

#' @rdname Reais
#' @examples
#' library(ggplot2)
#' p <- ggplot(centro_2015@data, aes(x = area_total, y = valor)) +
#'       geom_point(na.rm = TRUE)
#' p + scale_y_continuous(labels = reais(nsmall = 0))
#' @export
reais <- function(...) {
  function(x) Reais(x, ...)
}

#' Write in percentage form
#'
#' This function transforms any number to percentage format for reporting purposes.
#' @param x a real number
#' @examples
#' porcento(0.25)
#' pct(0.25)
#' @export
#'
porcento <- function (x, ...) {
  if (length(x) == 0)
    return(character())
  x <- brf(100*x, ...)
  paste0(x, "%")
}

#' @rdname porcento
#' @export
pct <- porcento

#' Regression equation
#'
#' Given a \link{lm} object, returns its regression equation
#' @param object object of class \code{lm}
#' @param type the equation type required: regression (reg) or estimation (est).
#' @param inline the equation mode. TRUE for inline equations or FALSE for
#' displayed mode.
#' @param FUN transformation applied to dependent variable.
#' @param accuracy number to round to; for POSIXct objects, a number of seconds
#' @param f rounding function: floor, ceiling or round
#' @examples
#' data(centro_2015)
#' centro_2015 <- within(centro_2015, PU <- valor/area_total)
#' fit <- lm(log(PU) ~ log(area_total) + quartos + suites + garagens +
#'            log(dist_b_mar) + padrao, centro_2015)
#' equacoes(fit)
#' equacoes(fit, accuracy = .01)
#' @export
equacoes <- function(object, type = c("reg", "est"), inline = TRUE, FUN,
                     accuracy = 100, f = round, errorTerm = TRUE){
  z <- object
  myformula <- stats::formula(z)
  parametros <- parameters(z)
  type <- match.arg(type)

  if(type == "reg") {
    lhs <- format(parametros$lhs)
  } else {
    lhs <- parametros$response
  }

  coefs <- coef(z)
  coefs <- plyr::round_any(coefs, accuracy, f)

  rhs <-  paste(coefs[1], "+", paste(coefs[-1], "\\cdot",
                                     names(coefs[-1]),
                                     collapse = " + ")
                )
  rhs <- gsub("\\_", "\\\\_", rhs)
  rhs <- gsub("\\+ -", "- ", rhs)
  rhs <- gsub("[.]", ",", rhs)

  if (type == "reg"){
    Formula <- paste(lhs, "=", rhs)
  } else if (!missing(FUN)) {
    Formula <- paste(lhs, " = ", appraiseR::invFunc(FUN), "(", rhs, ")",
                     sep = "")
  } else {
    message("Estimation regression asked but no transformation passed.")
  }

  if (type == "reg" & errorTerm == TRUE){
    Formula <- paste(Formula, "+ \\varepsilon")
  }

  if (inline == TRUE) {
    cat('$', Formula, '$', sep = "")
  } else {
    cat('$$', Formula, '$$', sep = "")
  }
}

#' Verify if a lm model has intercept
#'
#' Given a \link{lm} object, verifies if it has or not an intercept
#' @param object object of class \code{lm}
#' @examples
#' data(Alugueis)
#'
#' @export
hasIntercept <- function(obj, ...) attr(terms(obj, ...), "intercept")  == 1L

# adapted from `insight`
# jacobian / derivate for log models and other transformations ----------------
#' @export

aic <- function(object, adjust = TRUE) {
  if (adjust == FALSE) {
    x <- AIC(object)
  } else {
    x <- AIC(object)
    adjustment <- .ll_analytic_adjustment(object, insight::get_weights(object, na_rm = TRUE))
    x <- x - 2*adjustment
  }
  return(x)
}
#' @export
aicc <- function(object, adjust = TRUE) {
  k <- object$df.residual
  n <- length(object$residuals)
  p <- n - k - 1
  K <- 2 + p
  ll <- logLik(object)
  if (adjust == FALSE) {
    #x <- 2*K - 2*ll + 2*K*(K+1)/(n-K-1)
    x <- AIC(object) + 2*K*(K+1)/(n-K-1)
  } else {
    #x <- 2*K - 2*ll + 2*K*(K+1)/(n-K-1)
    x <- AIC(object) + 2*K*(K+1)/(n-K-1)
    adjustment <- .ll_analytic_adjustment(object, insight::get_weights(object, na_rm = TRUE))
    x <- x - 2*adjustment
  }
  return(x)
}

# this function adjusts any IC for models with transformed response variables
# (copied from insight)
#' @export
.adjust_ic_jacobian <- function(model, ic) {
  response_transform <- insight::find_transformation(model)
  if (!is.null(ic) && !is.null(response_transform) && !identical(response_transform, "identity")) {
    adjustment <- .safe(.ll_analytic_adjustment(model, insight::get_weights(model, na_rm = TRUE)))
    if (!is.null(adjustment)) {
      ic <- ic - 2 * adjustment
    }
  }
  ic
}


# adapted from `insight`
#' @export
.ll_analytic_adjustment <- function(x, model_weights = NULL) {
  tryCatch(
    {
      params <- parameters(x)
      trans <- params$depvarTrans
      resp <- na.omit(params$data[, params$response, drop = T])
      #resp <- insight::get_response(x)
      switch(trans,
             rec = .weighted_sum(log(100/resp^2), w = model_weights),
             rsqrt = .weighted_sum(log(5/resp^(3/2)), w = model_weights),
             identity = 0,
             log = .weighted_sum(log(1/resp), w = model_weights),
             log1p = .weighted_sum(log(1/(resp + 1)), w = model_weights),
             log2 = .weighted_sum(log(1/(resp * log(2))), w = model_weights),
             log10 = .weighted_sum(log(1/(resp * log(10))), w = model_weights),
             exp = .weighted_sum(resp, w = model_weights),
             expm1 = .weighted_sum((resp-1), w = model_weights),
             sqrt = .weighted_sum(log(0.5/ sqrt(resp)), w = model_weights),
             .ll_jacobian_adjustment(x, model_weights)
      )
    },
    error = function(e) {
      NULL
    }
  )
}


# this function calculates the adjustment for the log-likelihood of a model
# with transformed response (copied from insight)
#' @export
.ll_jacobian_adjustment <- function(model, weights = NULL) {
  tryCatch(
    {
      trans <- parameters(model)$depvarTrans
      .weighted_sum(log(
        diag(attr(with(
          insight::get_data(model, verbose = FALSE),
          stats::numericDeriv(
            expr = quote(trans(
              get(insight::find_response(model))
            )),
            theta = insight::find_response(model)
          )
        ), "gradient"))
      ), weights)
    },
    error = function(e) {
      NULL
    }
  )
}


.weighted_sum <- function(x, w = NULL, ...) {
  if (is.null(w)) {
    mean(x) * length(x)
  } else {
    stats::weighted.mean(x, w) * length(x)
  }
}