appraiseR: R Package For Real Estate Appraisals

#' Extracts grand or weighted grand means of a linear model with one or more
#' factors
#'
#' \code{grandMean} refits and object of class lm with new contrasts options
#' in order to obtain the grand mean or weighted grand mean of the response
#' variable when grouped by this (these) factor(s).
#' #' @param obj object of class lm
#'
#' @examples
#' # 1. Random generated data
#' # 1.1 One factor model
#' n = 25
#' set.seed(1)
#' Area <- rnorm(n = 3*n, mean = 360, sd = 50)
#' # Bairro <- gl(n = 3, k = n, labels = LETTERS[1:3])
#' Bairro <- sample(LETTERS[1:3], size = 3*n, replace = TRUE)
#'
#' dados <- data.frame(Area, Bairro = factor(Bairro))
#' dados <- within(dados, {
#'   PU <- ifelse(Bairro == "A", 7500,
#'                ifelse(Bairro == "B", 10000, 12500))
#'   PU <- PU - 5*Area + rnorm(n = 3*n, mean = 0, sd = 1000)
#' })
#'
#' summary(dados$Bairro) # unbalanced
#'
#' avg <- mean(dados$PU)
#' grandmean <- mean(tapply(dados$PU, dados$Bairro, mean))
#' weightedmean <- sum(tapply(dados$PU, dados$Bairro, mean)*summary(dados$Bairro))/75
#'
#' fit <- lm(PU ~ Bairro, data = dados)
#'
#' weightedMean(fit) # equals avg = weightedmean = mean(dados$PU)
#' grandMean(fit) # equals grandmean = mean(tapply(dados$PU, dados$Bairro, mean))
#'
#' all.equal(weightedMean(fit), mean(dados$PU)) # TRUE
#' all.equal(grandMean(fit), mean(tapply(dados$PU, dados$Bairro, mean))) # TRUE
#'
#' # 1.2 Same model, but whithout intercept
#'
#' # Things goes the same way for a model fitted without intercept
#'
#' fit <- lm(PU ~ Bairro - 1, data = dados)
#'
#' all.equal(weightedMean(fit), mean(dados$PU)) # TRUE
#' all.equal(grandMean(fit), mean(tapply(dados$PU, dados$Bairro, mean))) # TRUE
#'
#' # 1.3 Numeric covariate
#' # If we control for a numeric variable, the Weighted Mean will work the same
#'
#' fit1 <- lm(PU ~ Bairro + Area, data = dados)
#'
#' all.equal(weightedMean(fit1), mean(dados$PU)) # TRUE
#' all.equal(grandMean(fit1), mean(tapply(dados$PU, dados$Bairro, mean)))
#'
#' # 1.4 Model with two factors
#'
#' dados <- within(dados, AreaFactor <- cut(dados$Area, breaks = 3))
#'
#' fit2 <- lm(PU ~ Bairro + AreaFactor, data = dados)
#'
#' weightedMean(fit2)
#'
#' 2. Real data
#'
#' library(sf)
#' data(centro_2015)
#' centro_2015 <- within(centro_2015[1:50, ], PU <- valor/area_total)
#'
#' # Median centered covariates
#' fit <- lm(log(PU) ~ log(area_total/162) + I(quartos-3) + I(suites-1) +
#'           I(garagens-1) + log(dist_b_mar/595) + padrao, data = centro_2015)
#'
#' all.equal(grandMean(fit),
#'           mean(tapply(log(centro_2015$PU), centro_2015$padrao, mean))) #FALSE
#'
#' # Effects contrasts (show that intercept is the grand mean)
#'
#' fit1 <- lm(log(PU) ~ log(area_total/162) + I(quartos-3) + I(suites-1) +
#'           I(garagens-1) + log(dist_b_mar/595) + padrao, data = centro_2015,
#'           contrasts = list(padrao = "contr.sum"))
#'
#' Coefs <- fit1$coefficients
#' pbaixo <- Coefs['(Intercept)'] + Coefs['padrao1']
#' pmedio <- Coefs['(Intercept)'] + Coefs['padrao2']
#' palto <- Coefs['(Intercept)'] - Coefs['padrao1'] - Coefs['padrao2']
#' (pbaixo + pmedio + palto)/3 # equals the intercept of fit1
#'
#' # Weighted effects contrasts
#' library(wec)
#' contrasts(centro_2015$padrao) <- contr.wec(centro_2015$padrao, "baixo")
#'
#' fit2 <- lm(log(PU) ~ log(area_total/162) + I(quartos-3) + I(suites-1) +
#'           I(garagens-1) + log(dist_b_mar/595) + padrao, data = centro_2015)
#' summary(fit2)
#'
#' contrasts(centro_2015$padrao) <- contr.wec(centro_2015$padrao, "medio")
#'
#' fit3 <- lm(log(PU) ~ log(area_total/162) + I(quartos-3) + I(suites-1) +
#'           I(garagens-1) + log(dist_b_mar/595) + padrao, data = centro_2015)
#' summary(fit3)
#'
#' # Mean centered covariates
#' contrasts(centro_2015$padrao) <- contr.wec(centro_2015$padrao, "baixo")
#' fit4 <- lm(log(PU) ~ log(area_total/187.11) + I(quartos-2.66) + I(suites-1.2) +
#'           I(garagens-1.68) + log(dist_b_mar/545.52) + padrao,
#'           data = centro_2015)
#'
#' # Untouched covariates
#' contrasts(centro_2015$padrao) <- contr.wec(centro_2015$padrao, "baixo")
#' fit5 <- lm(log(PU) ~ log(area_total) + quartos + suites + garagens +
#'            log(dist_b_mar) + padrao, data = centro_2015)
#'
#' @export
weightedMean <- function(obj, ...){

  if (!hasIntercept(obj)) obj <- update(obj, ~.+1)

  residuos <- residuals(obj, "partial")
  k <- attr(residuos, "constant")
  k
}

#' Weighted Mean
#' @export
grandMean <- function(obj, ...){
  obj <- jtools::center_mod(obj)
  params <- parameters(obj)
  parametros <- params$parameters
  predictors <- params$predictors
  MF <- model.frame(obj)
  fatores <- parametros[which(sapply(MF, is.factor) == T)]

  contrs <- vector("list", length = length(fatores))
  names(contrs) <- fatores
  contrs[fatores] <- "contr.sum"

  if (hasIntercept(obj)){
    cl <- obj$call
    cl <- pryr::modify_call(cl, list(contrasts = contrs))
  } else {
    obj <- update(obj, ~.+1)
    cl <- obj$call
    cl <- pryr::modify_call(cl, list(contrasts = contrs))
  }

  newFit <- eval(cl, environment(stats::formula(obj)))
  WMean <- coef(newFit)["(Intercept)"]
  as.numeric(WMean)
}