Artigo.Rmd

---
title: "Crítica ao uso do Campo de Arbítrio do Avaliador"
subtitle:  "em situações de escassez de dados de mercado"
author: 
  - Luiz F. P. Droubi^[SPU/SC, lfpdroubi@gmail.com]
  - Carlos Augusto Zilli^[IFSC, carlos.zilli@ifsc.edu.br]
  - Willian Zonato^[SPU/SC, will.zonato@gmail.com]
  - Norberto Hochheim^[UFSC, hochheim@gmail.com]
date: "`r format(Sys.Date(), '%d/%m/%Y')`"
output:
  pdf_document:
    includes:
      in_header: preamble.tex
    keep_tex: yes
    latex_engine: xelatex
    number_sections: yes
    toc: no
  word_document: default
  html_document:
    fig_caption: yes
    keep_md: yes
classoption: a4paper, 11pt
documentclass: article
geometry: left=2.5cm,right=2.5cm,top=3cm,bottom=2.5cm
link-citations: yes
linkcolor: red
urlcolor: magenta
citecolor: green
csl: ABNT_UFPR_2011-Mendeley.csl
bibliography: skeleton.bib
always_allow_html: true
---

```{r, include = FALSE}
knitr::opts_chunk$set(echo = FALSE, fig.path = "images/", out.width = "65%",
                      dev = "CairoPNG", dpi = 600, fig.align = "center",
                      fig.pos = "H", warning = FALSE, message = FALSE)
type <- knitr::opts_knit$get("rmarkdown.pandoc.to")
library(appraiseR)
library(car)
library(knitr)
library(kableExtra)
library(xtable)
library(sjPlot)
library(plotly)
library(ggplot2)
library(cowplot)
theme_set(theme_cowplot(12))
library(dplyr)
library(lmtest)
#library(robustbase)
library(sandwich)
library(car)
library(mosaic)
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
red = gg_color_hue(1)
seed <- 1
epsilonsd <- 100
```

# Resumo {-}

Neste trabalho são apresentados aspectos teóricos e práticos relacionados ao
conceito de Campo de Arbítrio (CA) do Avaliador, dada a importância deste
conceito na Engenharia de Avaliações. Foram elencados os critérios previstos na
normativa que possibilitam ao avaliador fazer uso do Campo de Arbítrio,
detalhando cada um destes critérios levantados e ponderando se a adoção do
conceito de Campo de Arbítrio do avaliador é uma condição suficiente e
necessária para a solução dos problemas práticos enfrentados pelo avaliador.
Para melhor ilustrar, foram elaborados estudos de diversos casos com a geração
de dados randômicos simulando o problema da micronumerosidade de dados de uma
mesma característica, comparando os resultados obtidos com a adoção de diversas
abordagens, fazendo uso tanto do Campo de Arbítrio do Avaliador  quando do
intervalo de predição (IP) das previsões efetuadas com os modelos obtidos em
cada abordagem. Outro aspecto importante abordado lateralmente neste trabalho é
sobre a previsão de valores de venda a partir de dados de oferta, haja vista que
a falta de dados de transações e, em consequência a falta de um fator oferta
obtido cientificamamente, é um dos grandes motivos que levam os avaliadores a
fazerem uso do Campo de Arbítrio. Ao final, a partir da pesquisa elaborada e dos
resultados obtidos são feitas recomendações visando uma melhoria na NBR 14.653
numa eventual revisão desta.

# Introdução

No [XX Cobreap](http://www.cobreap.com.br/2019/), foi apresentado o artigo
intitulado "Crítica à Avaliação Intervalar na NBR14.653-02" [@droubi2019]. Neste
trabalho recomendou-se que, eventualmente, em uma próxima revisão da parte 2 da
NBR 14.653-02 [-@NBR1465302], a avaliação intervalar, tal como está normatizada
no momento, não seja mantida como se encontra. Com base nos resultados obtidos
naquele artigo, foram elaborados outros estudos, com ajustes diferentes, com o
intuito de verificar a pertinência das conclusões obtidas naquele trabalho,
assim como acrescentar alguns pontos no que tange ao Campo de Arbítrio do
avaliador, especialmente em relação à omissão de variável relevante. Também
foram feitas considerações a respeito do fator oferta. Pretende-se com este
trabalho deixar detalhes de uma proposta de como poder-se-ia modificar o atual
consenso da avaliação intervalar numa futura revisão normativa.

# Desenvolvimento e Fundamentação

Nesta seção serão abordados o conceito básico de campo de arbítrio, como
definido pela @NBR1465302, o conceito de viés devido à omissão de variável
relevante, tal como prevê o campo de arbítrio, além de questões correlatas, como
a questão do sobreajuste (*overfitting*) em modelos estatísticos e o *tradeoff*
entre viés e variância.

## Campo de Arbítrio do Avaliador

A NBR 14.653-02 [@NBR1465302] define que:

> "_o Campo de Arbítrio pode ser utilizado quando **variáveis relevantes** para
a avaliação do imóvel **não tiverem sido contempladas no modelo**, por escassez
de dados de mercado, por inexistência de fatores de homogeneização aplicáveis ou
porque essas variáveis não se apresentaram estatisticamente significantes em
modelos de regressão, **desde que a amplitude de até mais ou menos 15% seja
suficiente para absorver as influências não consideradas e que os ajustes sejam
justificados**_".

Um ponto levantado por Droubi *et al.* [-@droubi2019] é que a comparação dos
valores arbitrados dentro do CA não deveriam ser comparados ao intervalo de
confiança da regressão linear, pois *o intervalo de confiança é para a média*,
ou seja, o intervalo de confiança é apenas para aferir a precisão do cálculo dos
valores ajustados pela reta de regressão (média + variância explicada) e não
serve para contemplar outros efeitos não calculados pelo modelo para a avaliação
do bem-avaliando. Para isto, a comparação devida seria com o intervalo de
predição, ou seja, o intervalo onde se encontram os valores reais observados do
mercado e não apenas os valores médios (a soma da média e da variância explicada
e não-explicada).

Este primeiro aspecto refere-se mais à questão da avaliação intervalar do que ao
campo de arbítrio em si, e poderia ser contornado simplesmente com a mudança dos
critérios da avaliação intervalar, adotando-se critérios mais adequados para a 
elaboração de intervalos.

Um outro ponto igualmente importante levanto por Droubi *et al.* [-@droubi2019],
porém, diz respeito apenas ao CA: o intervalo fixo, de $\pm$ 15%, pode até ser
razoável em algumas situações, mas pode levar a valores de baixíssima
probabilidade de ocorrência prática em mercados de baixa variabilidade, ou
ainda, pelo contrário, num mercado de grande variância, levar a valores que
sequer chegam próximos dos valores extremos observados no mercado.

Porém, entende-se que Droubi *et al.* [-@droubi2019] não esclareceram
completamente os motivos pelos quais um modelo pode vir a apresentar um
intervalo de predição tão grande, enquanto noutros o intervalo de predição pode
ser mais estreito.

Como será visto a seguir, isto pode ser uma questão relacionada apenas ao
próprio mercado, como defenderam Droubi *et al.* [-@droubi2019], mas também pode
estar relacionado às variáveis importantes não contempladas no modelo. Na
próxima subseção será analisado um dos fatores que podem levar à omissão de
variáveis importantes num modelo de avaliação de imóveis.

## Escassez de dados de mercado - Micronumerosidade

Um dos fatores que levam à omissão de variáveis relevantes do modelo, o que 
embasa a utilização do CA do avaliador, é o critério da micronumerosidade.

De acordo com o Anexo A da referida norma, quando da utilização de variáveis
dicotômicas ou qualitativas, um número mínimo de dados de cada característica,
variando de 3 a 10, a depender do tamanho da amostra, devem ser efetivamente
utilizados. Em geral, o número de dados de cada característica ($n_i$) deve ser
igual a 10% do número total de dados ($n_i \geq 10\% \ n$) da amostra, limitado
a um mínimo de 3 dados caso a amostra possua um número total de dados menor do
que 30 ($n<30$) e um limitante superior de 10 dados $n_i \geq 10$ caso a amostra
possua um número de dados maior do que 100 ($n > 100$).

A NBR 14.653-02 não é clara, no entanto, sobre como proceder no caso de o
avaliador encontrar em sua amostra um número de dados de uma característica
menor do que o mínimo estabelecido. Existem diversas possibilidades, como:

a. a remoção da variável que apresenta a micronumerosidade;
b. a recodificação da variável (quando possível);
c. a remoção dos dados que apresentam aquela característica em número 
insuficiente. 

A depender do procedimento adotado, o avaliador obterá um modelo diferente e a
implicância da adoção de cada procedimento será mostrada nos estudos de casos
apresentados. Presume-se, no entanto, a partir de uma leitura global da norma,
que a mesma sugere que a variável **e** os dados sejam removidos do modelo, pois
com a remoção apenas da variável e a inclusão dos dados com características
diferentes no modelo, a hipótese da homoscedasticidade dos resíduos,
normalmente, não se verifica, como será mostrado, e a norma também veda a
utilização de modelos em que as hipóteses da inferência clássica não se cumpram,
de acordo com o *caput* do item A.2.

## Inexistência de fatores de homogeneização

Ficou claro durante o [XX Cobreap](http://www.cobreap.com.br/2019/) que grande
parte dos avaliadores se utiliza do CA para resolver um dos maiores problemas da
prática atual da Engenharia de Avaliações, que é a falta de dados de transações,
o que implica na impossibilidade do cálculo de um fator oferta. A imensa maioria
dos Avaliadores se valem de uma pesquisa de dados ofertados, chegando assim a um
modelo para os preços de oferta, que infelizmente não corresponde a um modelo de
preços de venda dos imóveis.

A atual normativa não permite a aplicação de um fator oferta aos dados
pesquisados, a menos que este fator seja deduzido com a utilização de
metodologia científica.

No entanto, um modelo econométrico proposto por Horowitz [-@horowitz] permite a
previsão de preços de venda a partir de preços de oferta mais precisamente do
que os que seriam obtidos com um modelo hedônico elaborado à partir dos próprios
dados de venda.

Num estudo de caso Horowitz [-@horowitz, 124-125], encontrou um modelo
econométrico que previu os preços de vendas a partir dos preços de oferta com um
erro médio quadrático de apenas \$2.960, enquanto que com um modelo hedônico
elaborado a partir dos próprios dados de venda foi obtido um erro de \$10.361,
ou seja, um erro médio muito maior do que o erro médio obtido com o modelo
econométrico proposto.

Está fora do escopo deste trabalho detalhar a abordagem de Horowitz
[-@horowitz]. Contudo, a existência de um método capaz de prever preços de venda
de imóveis a partir apenas de dados de oferta, e com melhor precisão do que um
modelo hedônico elaborado com os próprios dados de venda é um fator a menos em
favor da utilização do CA do avaliador com esta finalidade. Estudos devem ser
levados a cabo e possivelmente o modelo tenha que ser adaptado para se ajustar à
realidade da prática da Engenharia de Avaliações no Brasil, porém a existência
de tal abordagem não pode permanecer completamente ignorada como é hoje na
Engenharia de Avaliações.

## Falta de significância dos regressores

Um dos pontos polêmicos da NBR 14.653-02 consiste na adoção de graus de 
fundamentação de modelos de regressão linear. Entre os critérios utilizados para
o enquadramento dos modelos nos diversos graus de fundamentação está o nível de
significância máxima para a rejeição da hipótese nula de cada regressor. Está
além do escopo deste artigo uma crítica pormenorizada a este critério, porém
deve-se salientar que este tipo de análise de importância de regressores à
partir dos p-valores dos testes de hipótese tem sido muito contestada entre os
estatísticos [ver @matloff2017, 349; @ASA].

## Viés devido à omissão de variável

Nesta subseção serão apresentados aspectos teóricos a respeito da omissão de 
variáveis relevantes de um modelo de regressão, independentemente dos motivos que
levem a esta omissão, seja por escassez de dados de mercado, seja por falta de
significância do regressor.

Segundo Matloff [-@matloff2017, 24], adicionando mais variáveis independentes a
um modelo se está a reduzir o viés do modelo.

Em tese, qualquer coisa que influencie o valor da variável dependente e não
esteja presente no modelo é absorvido pelo termo de erro (variação não-explicada
pelo modelo), que se assume estar não-correlacionado com os regressores
($\mathbb{E}[\varepsilon|X] = 0$). No entanto, se houver algum grau de correlação
entre a variável omitida e as variáveis presentes, como os efeitos da variável
omitida são absorvidas pelo termo de erro, a hipótese de independência dos 
resíduos não será mantida. Segundo Matloff [-@matloff2017, 25], um modelo assim
é chamado pelos economistas de modelo mal-especificado (*misspecified*), 
enquanto os estatísticos tratam este de tipo viés como viés de modelagem (*model bias*).

De acordo com Matloff [-@matloff2017, 25], no entanto, isto não quer dizer que
se trata de um viés que tenha sido adicionado deliberadamente ao modelo, apenas
quer dizer que existe no modelo um viés sistêmico, que é permanente e não
desaparece, por mais que aumente-se o tamanho da amostra.

## Sobreajuste 

O sobreajuste ou *overfitting* em estatística é uma prática que consiste em
elaborar um modelo qualquer de maneira tão complexa que, ao invés de captar o
sinal, capta o ruído [@matloff2017, 24]. 

A prática, então, de se adicionar variáveis ao modelo para reduzir o viés da
estimação, como visto na subseção anterior, tem um limite: a introdução de uma
variável, ao diminuir o viés, também aumenta a variância das estimativas dos
coeficientes do modelo. Em um determinado ponto, o benefício da adição de 
variáveis ao modelo será menor do que a maior variabilidade que ocorrerá na 
estimação, o que é chamado de *overfitting* [@matloff2017, 24].

Segundo Matloff [-@matloff2017, 24], para qualquer estimador estatístico $\hat
\theta$ com variância finita, pode ser mostrado que:

\begin{equation} \label{eq:tradeoff}
MSE(\hat \theta) = Var(\hat \theta) + B^2(\hat \theta)
\end{equation}

Onde $MSE(\hat\theta)$ é o erro médio quadrático associado ao estimador, 
$Var(\hat\theta)$ é a variância do estimador e $B^2(\hat\theta)$ é o quadrado do
viés do estimador.

Segundo Matloff [-@matloff2017, 25], a equação \ref{eq:tradeoff} mostra que,
adicionando variáveis ao estimador de regressão linear, reduz-se o termo
$B^2(\hat \theta)$, porém ao custo do aumento da variância ($Var(\hat\theta)$) e
vice-versa, ao se remover variáveis do modelo, se está aumentando o viés do 
modelo, reduzindo, porém sua variância. De acordo com Matloff pode ser até
benéfico remover variáveis de um modelo com o intuito de reduzir a variância do
estimador, ao custo de aceitar um pouco de viés [@matloff2017, 25].

Um bom exemplo para o entendimento do problema pode ser encontrado em
@matloff2017 [p. 342-343], reproduzido abaixo:

Supondo uma amostra com alturas de meninos e meninas, estimar as alturas médias
de ambas as populações, com variância conhecida igual a $\sigma^2$. De acordo 
com Matloff, os estimadores naturais seriam:

\begin{equation} \label{eq:boysandgirls}
\hat\mu_1 = \frac{1}{n}\sum_{i=1}^{n}X_i \qquad e \qquad \hat\mu_2 = \frac{1}{n}\sum_{i=1}^{n}Y_i
\end{equation}

Porém, se o número de dados da amostra é pequeno, pode ser melhor simplificar o 
modelo e utilizar um único estimador para ambos os grupos:

\begin{equation} \label{eq:all}
\check \mu_i = \frac{1}{2}(\bar X + \bar Y), \ i = 1,2
\end{equation}

O critério para a escolha do melhor estimador é o erro médio quadrático (MSE).

Ocorre que o erro médio quadrático para os estimadores das equações 
\ref{eq:boysandgirls} é igual à variância dos estimadores, dado que 
ambos são estimadores não-viesados ($B(\hat \mu_1) = B(\hat \mu_2) = 0$):

\begin{equation} \label{eq:mse1}
MSE(\hat\mu_1) + MSE(\hat\mu_2) = 2 Var(\hat \mu_i) = 2 \frac{\sigma^2}{n}
\end{equation}

Já para os estimadores $\check \mu_i$ a soma dos seus erros médios quadráticos
pode ser vista na equação \ref{eq:mse2} [^1]

\begin{equation} \label{eq:mse2}
MSE(\check\mu_1) + MSE(\check\mu_2) = \frac{\sigma^2}{n} + \frac{1}{2}(\mu_1 - \mu_2)^2
\end{equation}

[^1]: Obtido através da soma das variâncias e viés quadrados dos estimadores, 
como pode ser visto em @matloff2017 [p. 386-387].

Desta maneira, para saber qual o melhor estimador, deve-se comparar os valores
dos erros médios quadráticos obtidos nas equações \ref{eq:mse1} e \ref{eq:mse2}.

Caso $2\frac{\sigma^2}{n} < (\mu_1-\mu_2)^2$, adotam-se os estimadores naturais. 
Caso contrário, mais vale a pena adotar o estimador único para ambos os grupos.

Deve-se ter em conta que:

1. se a diferença entre as médias $\mu_1$ e $\mu_2$ é pequena, o estimador único é bem próximo da realidade;
2. se o número de dados $n$ é pequeno, não há dados suficientes para se estimar com precisão os dois estimadores em separado;
3. se a variância da população é alta, também será preciso mais dados para a estimação em separado das duas médias.

Segundo Matloff [-@matloff2017, p. 344], este exemplo é um problema trivial de
regressão, onde deve-se escolher entre um modelo nulo ($p = 0$), e um modelo 
mais elaborado, com um regressor ($p = 1$), a saber, a variável dicotômica que 
indica se tratar de um dado de altura de menino ou menina.

Ainda de acordo com Matloff [-@matloff2017, p. 344], isto pode ser estendido ao
contexto genérico de modelos de regressão e classificação, com $p$ regressores e
$n$ observações: se $p$ é grande e/ou $n$ é pequeno, pode ser desejável utilizar
um modelo mais simples, omitindo alguns regressores em que o valor de $\beta$ 
seja negligenciável.

De acordo com Tukey e outro trabalho mais recente de Portnoy [*apud*
@matloff2017, p.344], via de regra, devem existir no máximo $\sqrt{n}$
regressores em um modelo de regressão para se evitar o sobreajuste do modelo.

# Estudos de Casos

Neste trabalho foram elaborados dois estudos de casos. 

No primeiro estudo de caso, que visa explorar os efeitos da omissão de variáveis
relevantes nos modelos de regressão linear, foram feitas quatro simulações,
variando o efeito da majoração dos valores de lotes em situação de esquina, em
relação aos lotes em situação de meio de quadra.

Foi elaborado ainda um segundo estudo de caso visando obter reflexões acerca da
utilização do CA quando, pelo contrário, não existem variáveis **relevantes** 
omitidas, ainda que haja variabilidade na amostra, o que também é normal que 
aconteça.

## Estudo de Caso 1 -- Viés devido à omissão de variável relevante

Neste Estudo de Caso foram elaborados quatro simulações de dados aleatórios para 
discussão de alguns fatos que são ilustrados com a aplicação de diferentes 
abordagens para o tratamento dos dados.

Para os quatro casos, os dados do Valor Unitário de cada lote (VU) foram 
criados através da seguinte equação de regressão (da população), conforme
equação \ref{eq:generic}:

\begin{equation} \label{eq:generic}
VU = \beta_0 + \beta_1 \cdot Area + \beta_2 \cdot Situacao + \varepsilon
\end{equation}

Onde $\varepsilon \sim \mathcal{N}(0, \sigma^2)$. Para os diversos casos foram
variados os valores dos coeficientes da população $\beta_i$ utilizados para a
geração de dados, de maneira a salientar o comportamento dos diversos modelos
ajustados tanto na estimação dos coeficientes ($\hat\beta_i$), quando na 
elaboração de previsões sobre um lote de referência.

A vantagem de se trabalhar com dados gerados é que se conhece *a priori* o valor
dos coeficientes ($\beta_i$) da população, já que estes estão pré-definidos.
Obviamente que, com a inclusão de um ruído ($\varepsilon$), no ajuste dos modelos
de regressão aos dados gerados, não são obtidos para os coeficientes ajustados
($\hat \beta_i$) exatamente os mesmos valores dos coeficientes populacionais
($\beta_i$) utilizados para a geração dos dados. Mas se o modelo estiver bem
especificado, os valores estimados para os coeficientes devem ser próximos dos
coeficientes utilizados para a geração dos dados, ou seja, $\hat \beta_i \approx
\beta_i$.

Em todos os casos foram gerados dados (propositalmente) em número insuficiente 
($<3$) para o caso dos lotes em situação de esquina.

Foram estudadas, então, as consequências da utilização de três diferentes tipos
de abordagens:

- O ajuste de um modelo de regressão com a inclusão de todos os dados, porém com
a omissão da variável `situacao`, visando emular o caso previsto na definição de 
CA, de omissão de variável relevante do modelo devido à escassez
de dados de mercado;

- O ajuste de um modelo de regressão com a exclusão dos dados de esquina, de 
maneira a se obter um modelo ajustado em cima de uma amostra homogênea, isto é,
um amostra onde todos os dados possuem as mesmas características em relação à
variável `situacao`.

- O ajuste de um modelo de regressão com a inclusão de todos os dados e todas as
variáveis, apesar da recomendação em contrário da NBR 14.653-02.

### Mercado 1

```{r dadosEx1, echo = FALSE}
set.seed(seed)
area <- c(360, 360, 360, 360, 480, 480, 360, 360, 480, 480)
situacao <- c(rep(0, 4), rep(1, 2), rep(0, 4))
VU <- 5000 - 5*area + 250*situacao + rnorm(10, 0, epsilonsd)
dados <- data.frame(VU, Area = area, 
                    Situacao = factor(situacao, 
                                      levels = c(0, 1), 
                                      labels = c("Meio de Quadra", "Esquina")
                                      )
)
```

```{r}
xtable(dados, 
       caption = "Dados para o M. 1.",
       label = "tab:ex1",
       digits = 0,
       align = rep("c", 4)) %>%
  xtable2kable(include.rownames = FALSE) %>%
  kable_styling(latex_options = "striped",
                position = "float_right")
```

Para este exemplo os dados foram gerados através da equação \ref{eq:ex1}, onde
$\varepsilon \sim \mathcal{N}(0, \,100^2)$:

\begin{equation} \label{eq:ex1}
VU = 5000 - 5 \cdot Area + 250 \cdot Situacao + \varepsilon
\end{equation}

Seja o caso de se avaliar um lote urbano de esquina com 480$m^2$, baseado na
amostra exibida na tabela \ref{tab:ex1}, obtida com as simulações, onde
`Situacao` é uma variável dicotômica que diferencia os imóveis em situação de
meio de quadra (0) dos imóveis de esquina (1).

De acordo com a NBR 14.653-02, como estão disponíveis apenas 2 dados de mercado
em situação de esquina, não seria possível a utilização da variável `Situacao`,
devido à micronumerosidade. Neste tipo de situação (escassez de dados de
mercado), a norma permite a utilização do CA.

Ou seja, segundo a NBR 14.653-02 leva a crer, poderia ser elaborado um modelo
com todos os dados amostrais, com a inclusão apenas da variável `Area`,
omitindo-se a variável `Situacao` e utilizar o CA para majorar o valor estimado
do lote pelo modelo, por este estar em situação de esquina, *desde que a
amplitude de até mais ou menos 15% seja suficiente para absorver as influências
não consideradas e que os ajustes sejam justificados*.

```{r, echo = FALSE}
#fit_a <- lmrob(VU ~ Area, data = dados)
fit_a <- lm(VU ~ Area, data = dados)
fit1_a <- lm(VU ~ Area, dados, subset = -c(5,6), x = TRUE, y = TRUE)
fit2_a <- lm(VU ~ Area + Situacao, dados, x = TRUE, y = TRUE)
```

```{r}
# Adjust standard errors
cov1         <- vcovHC(fit_a, type = "HC3")
robust_se    <- sqrt(diag(cov1))
```

```{r, results='asis', echo = FALSE}
stargazer::stargazer(fit_a, fit1_a, fit2_a, label = "tab:tab1", 
                     title = "Comparação de modelos para o Ex. 1.",
                     type = type, header = FALSE, table.placement = "H",
                     ci = TRUE, ci.level = .80, report = "vcst*", single.row = T,
                     intercept.bottom = FALSE, intercept.top = TRUE,
                     decimal.mark = ",", digit.separator = ".",
                     digits = 1, star.cutoffs = c(0.30, 0.20, 0.10),
                     # se = list(robust_se, NULL, NULL),
                     # notes = c("Reportados erros robustos para modelo da coluna 1."),
                     notes.label = "Notas:", font.size = "footnotesize")
# tab_model(fit_a, fit1_a, fit2_a,
#           dv.labels = c("Modelo 1", "Modelo 2", "Modelo 3"),
#           show.ci = 0.80,
#           show.stat = TRUE,
#           show.intercept = TRUE)
```

Na tabela \ref{tab:tab1} são mostrados os coeficientes dos três modelos
ajustados, isso fica claro: no primeiro modelo (coluna 1), foram utilizados
todos os dados e apenas a variável `Area`; no segundo modelo (coluna 2), foram
utilizados apenas os dados em situação de meio-de-quadra; e no terceiro modelo
(coluna 3) foram utilizados todos os dados e foi incluída a variável `Situacao`,
apesar da micronumerosidade.

O que se percebe é que o primeiro modelo viesou para cima o coeficiente da 
variável `Area`. Ou seja, o efeito da `Situacao` do lote, que não foi incluso
no modelo, foi "absorvido" pela outra variável restante, ou seja, a variável 
`Area`. Isto ocorre porque a modelagem não respeitou um dos princípios da 
inferência clássica, que é a homogeneidade da amostra, haja vista que a 
amostra utilizada mistura imóveis em situação de meio de quadra e imóveis em
situação de esquina, sem uma variável que modele esta diferença entre os dados
amostrais.

Os outros dois modelos se aproximaram melhor dos coeficientes "reais" de 
regressão, ou seja, o valor dos coeficientes para a população ( 5000, -5,0 e
250). Isto se deve à homogeneidade da amostra, obtida ou através da exclusão dos
dados não-homogêneos, ou com a inclusão de todos os dados e variáveis
relevantes, em que pese estar em desacordo com o estabelecido pela NBR 14653-02.
Deve-se notar, contudo, que a estimação foi melhor para os coeficientes obtidos
com o terceiro modelo, haja vista o maior número de dados.

Pelo lado da estimação  dos coeficientes, portanto, não se recomenda a
utilização da primeira abordagem. A segunda abordagem é preferível à primeira,
porém a terceira abordagem é ainda mais precisa que a segunda.

Com a primeira abordagem (o modelo da coluna 1), onde foram misturados dados de
esquina e de meio de quadra, esperava-se desvios em relação à homoscedasticidade
dos resíduos. No entanto, provavelmente devido à pequena diferença entre os
valores dos lotes em situação de esquina e dos lotes em situação de meio de
quadra, a hipótese da homoscedasticidade não pôde ser rejeitada via teste de
Breusch-Pagan:

```{r}
bptest(fit_a)
```

```{r}
newx <- seq(min(dados$Area), max(dados$Area), by=2.5)
pred.int1 <- predict(fit_a, newdata = data.frame(Area = newx), 
                         interval = "prediction", level = 0.80)
pred.int1 <- data.frame(VU = pred.int1, Area = newx)
```

```{r, echo = FALSE}
# p_a <- Predict(fit_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
#                     interval = "confidence", level = .80, vcov.= hccm(fit_a))
p_a <- predict(fit_a, newdata = data.frame(Area = 480, Situacao = "Esquina"),
                    interval = "confidence", level = .80)
p_a_1 <- predict(fit_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
                 interval = "confidence", level = 0.80)
p1_a <- predict(fit1_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
                interval = "confidence", level = 0.80)
p2_a <- predict(fit2_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
                interval = "confidence", level = 0.80)
# P_a <- Predict(fit_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
#                interval = "prediction", level = 0.80, vcov.= hccm(fit_a))
P_a <- Predict(fit_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
               interval = "prediction", level = 0.80)
P1_a <- predict(fit1_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
                interval = "prediction", level = 0.80)
P2_a <- predict(fit2_a, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
                interval = "prediction", level = 0.80)
amp_a <- amplitude(P_a)
amp1_a <- amplitude(P1_a)
amp2_a <- amplitude(P2_a)
ca_a <- campo_arbitrio(p_a)
ca1_a <- campo_arbitrio(p1_a)
ca2_a <- campo_arbitrio(p2_a)
```

```{r figex1, fig.cap="Gráfico dos modelos para o primeiro mercado."}
plt1 <- ggplot(pred.int1, aes(x = Area, y = VU.fit)) +
  geom_point(data = dados, aes(x = Area, y = VU, color = Situacao), size = .5) +
  geom_smooth(data = pred.int1, aes(ymin = VU.lwr, ymax = VU.upr), 
              stat = "identity", size = .5, color = "orange") +
  geom_abline(intercept = coef(fit1_a)[1], slope = coef(fit1_a)[2], 
              color = "blue", size = .5) +
  geom_segment(aes(x = 360, y = 1.15*pred.int1[1, 1], 
                   xend = 480, yend = 1.15*pred.int1[49, 1]
                   ),
               color = "orange", lty = "ff", size = .5) +
  geom_segment(aes(x = 360, y = 0.85*pred.int1[1, 1], 
                   xend = 480, yend = 0.85*pred.int1[49, 1]
                   ),
               color = "orange", lty = "ff", lwd = .5) +
  scale_y_continuous(name = "VU", expand = c(0.025, 0.025), limits = c(2000, 4500)) +
  scale_x_continuous(expand = c(0.025, 0.025), limits = c(360, 480)) +
  scale_color_hue(name = 'Situação', labels = c('Meio de Quadra', 'Esquina')) +
  theme(legend.position = 'bottom', legend.direction = "horizontal",
        axis.text.x = element_text(color = "grey20", size = 8),
        axis.text.y = element_text(color = "grey20", size = 8))
```

A tabela 3 apresenta o comportamento das diversas abordagens para a previsão de
novos valores com a utilização dos modelos ajustados. Nesta tabela podem ser
vistos os valores para a estimativa central e para os limites inferiores e
superiores do IP (\@80%) e do CA, assim como os limites do CA. O valor teórico, 
obtido utilizando-se os coeficientes de regressão para a população ($\beta_i$) 
seria de `r reais()(5000 - 5*480 + 250)`.

Percebe-se desta maneira que para o **primeiro modelo**, o valor adicional do
lote em esquina já foi absorvido e já se encontra, em grande parte, na
estimativa central efetuada com o mesmo. No **segundo modelo**, que é um modelo
feito apenas para os dados de lotes em situação de meio de quadra, não há
qualquer influência de lotes em situação de esquina, ficando a previsão de
valores para lotes de esquina comprometida. E o **terceiro modelo**, apesar da
micronumerosidade, foi o que melhor estimou o lote em situação de esquina.

| Abordagem | Estimativa central     |$IP_{inf}$             | $IP_{sup}$              | Amplitude IP (%)    | $CA_{inf}$       | $CA_{sup}$        |
|:----------|-----------------------:|----------------------:|------------------------:|----------------:|-----------------:|------------------:|
| 1         |`r brf(p_a[, "fit"])`   |`r brf(P_a[, "lwr"])`  |  `r brf(P_a[, "upr"])`  | `r brf(amp_a)`  | `r brf(ca_a[1])` | `r brf(ca_a[2])`  |  
| 2         |`r brf(p1_a[, "fit"])`  | `r brf(P1_a[, "lwr"])`|  `r brf(P1_a[, "upr"])` | `r brf(amp1_a)` | `r brf(ca1_a[1])`| `r brf(ca1_a[2])` |
| 3         |`r brf(p2_a[, "fit"])`  | `r brf(P2_a[, "lwr"])`|  `r brf(P2_a[, "upr"])` | `r brf(amp2_a)` | `r brf(ca2_a[1])`| `r brf(ca2_a[2])` |
Table: Comparação dos limites do IP e CA para os modelos do Ex. 1.

Deve-se notar que, para este caso, o CA é *suficiente para 
para absorver as influências não consideradas*, ou seja, com a adoção do limite
superior do CA é possível chegar e ultrapassar o valor 'real' do
imóvel, de R\$2.800,00.

### Mercado 2

```{r dadosEx2, echo = FALSE}
set.seed(seed)
area <- c(360, 360, 360, 360, 480, 480, 360, 360, 480, 480)
situacao <- c(rep(0, 4), rep(1, 2), rep(0, 4))
VU <- 5000 - 5*area + 390*situacao + rnorm(10, 0, epsilonsd)
dados <- data.frame(VU, Area = area, 
                    Situacao = factor(situacao, 
                                      levels = c(0, 1), 
                                      labels = c("Meio de Quadra", "Esquina"))
)
```

```{r}
xtable(dados, 
       caption = "Dados para o M. 2.",
       label = "tab:ex2",
       digits = 0,
       align = rep("c", 4)) %>%
  xtable2kable(include.rownames = FALSE) %>%
  kable_styling(latex_options = "striped", 
                full_width = F, 
                position = "float_right")
```

Neste segundo exemplo, foram seguidos exatamente os mesmos passos do exemplo 
anterior. Contudo, para a geração dos valores, foi utilizado um peso maior para
os lotes de esquina, conforme a equação \ref{eq:ex2}, onde 
$\varepsilon \sim \mathcal{N}(0, \,100^2)$ :

\begin{equation} \label{eq:ex2}
VU = 5000 - 5 \cdot area + 390 \cdot situacao + \varepsilon
\end{equation}

Os dados gerados encontram-se na tabela \ref{tab:ex2}. Deve-se reparar que,
aqui, erros à parte, foi adicionado R$390,00 ao valor de um lote em situação de
esquina em relação a um lote em situação de meio de quadra, o que significa um
adicional de exatamente 15% ao valor do lote-padrão de 480 $m^2$ em situação de
meio de quadra, que apresenta valor 'real' de R\$2.600,00. Ou seja, o valor do
lote que se pretende prever, de 480 $m^2$, em situação de esquina, neste mercado, 
terá valor 'real' de `r reais()(5000 - 5*480 + 390)`.

```{r, echo = FALSE}
# fit_b <- lmrob(VU ~ Area, data = dados)
fit_b <- lm(VU ~ Area, data = dados)
fit1_b <- lm(VU ~ Area, dados, subset = -c(5,6), x = TRUE, y = TRUE)
fit2_b <- lm(VU ~ Area + Situacao, dados, x = TRUE, y = TRUE)
```

Novamente, então, foram ajustados três modelos, como os descritos no exemplo
anterior.

Na tabela \ref{tab:tab2} são mostrados os coeficientes dos três modelos
ajustados. Novamente fica claro que: no primeiro modelo (coluna 1), os 
coeficientes ajustados estão viesados, no segundo e no terceiro modelos eles 
estão estimados corretamente, porém o terceiro modelo é mais preciso (menores
intervalos de confiança).

```{r}
# Adjust standard errors
cov1         <- vcovHC(fit_b, type = "HC3")
robust_se    <- sqrt(diag(cov1))
```

```{r, results='asis', echo = FALSE}
stargazer::stargazer(fit_b, fit1_b, fit2_b, label = "tab:tab2",
                     title = "Comparação de modelos para o Ex. 2.",
                     type = type, header = FALSE, table.placement = "H",
                     ci = TRUE, ci.level = .80, report = "vcst*", single.row = T,
                     intercept.bottom = FALSE, intercept.top = TRUE,
                     decimal.mark = ",", digit.separator = ".",
                     digits = 1, star.cutoffs = c(0.30, 0.20, 0.10),
                     se = list(robust_se, NULL, NULL),
                     notes = c("Reportados erros robustos para modelo da coluna 1."),
                     notes.label = "Notas:", font.size = "footnotesize")
# tab_model(fit_b, fit1_b, fit2_b,
#           dv.labels = c("Modelo 1", "Modelo 2", "Modelo 3"),
#           show.ci = 0.80,
#           show.stat = TRUE,
#           show.intercept = TRUE)
```

```{r, echo = FALSE}
p_b <- Predict(fit_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80, vcov. = hccm(fit_b))
#p_b <- predict(fit_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
 #              interval = "confidence", level = 0.80)
p1_b <- predict(fit1_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
p2_b <- predict(fit2_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
P_b <- Predict(fit_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80, vcov. = hccm(fit_b))
#P_b <- predict(fit_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
#               interval = "prediction", level = 0.80)
P1_b <- predict(fit1_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
P2_b <- predict(fit2_b, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
amp_b <- amplitude(P_b)
amp1_b <- amplitude(P1_b)
amp2_b <- amplitude(P2_b)
ca_b <- campo_arbitrio(p_b)
ca1_b <- campo_arbitrio(p1_b)
ca2_b <- campo_arbitrio(p2_b)
```

Com o aumento da diferença de valores médios para os lotes em situação de
esquina em relação aos lotes em situação de meio de quadra, com o modelo da
primeira abordagem a homoscedasticidade não pode ser verificada, segundo o teste
de Breusch-Pagan:

```{r}
bptest(fit_b)
```

O intervalo de predição \@80% e os limites do CA para os três modelos na tabela
6 [^2]. Deve-se notar que, em relação ao primeiro estudo de caso, a amplitude do IP
para o primeiro modelo aumentou, enquanto a amplitude do IP para o terceiro
modelo diminuiu. Na segunda abordagem, o que se obtêm é o mesmo modelo do
primeiro estudo de caso, portanto os intervalos de predição e de CA são
exatamente os mesmos do caso anterior.

[^2]: O intervalos de predição para o modelo da primeira abordagem foram 
calculados com a utilização de erros robustos.

```{r}
newx <- seq(min(dados$Area), max(dados$Area), by=2.5)
pred.int2 <- predict(fit_b, newdata = data.frame(Area = newx), 
                         interval = "prediction", level = 0.80)
pred.int2 <- data.frame(VU = pred.int2, Area = newx)
```

```{r figex2, fig.cap="Gráfico dos modelos para o segundo mercado."}
plt2 <- ggplot(pred.int2, aes(x = Area, y = VU.fit)) +
  geom_point(data = dados, aes(x = Area, y = VU, color = Situacao), size = .5) +
  geom_smooth(data = pred.int2, aes(ymin = VU.lwr, ymax = VU.upr), 
              stat = "identity", size = .5, color = "orange") +
  geom_abline(intercept = coef(fit1_b)[1], slope = coef(fit1_b)[2], 
              color = "blue", size = .5, lwd = 1) +
  geom_segment(aes(x = 360, y = 1.15*pred.int2[1, 1], 
                   xend = 480, yend = 1.15*pred.int2[49, 1]
                   ),
               color = "orange", lty = "ff") +
  geom_segment(aes(x = 360, y = 0.85*pred.int2[1, 1], 
                   xend = 480, yend = 0.85*pred.int2[49, 1]
                   ),
               color = "orange", lty = "ff") +
  scale_y_continuous(name = "VU", expand = c(0.025, 0.025), limits = c(2000, 4500)) +
  scale_x_continuous(expand = c(0.025, 0.025), limits = c(360, 480)) +
  theme(legend.position = 'bottom', legend.direction = "horizontal",
        axis.text.x = element_text(color = "grey20", size = 8),
        axis.text.y = element_text(color = "grey20", size = 8))
```


| Modelo | Estimativa central     |$IP_{inf}$             | $IP_{sup}$              | Amplitude IP (%)    | $CA_{inf}$       | $CA_{sup}$        |
|:-------|-----------------------:|----------------------:|------------------------:|----------------:|-----------------:|------------------:|
| 1      |`r brf(p_b[, "fit"])`   |`r brf(P_b[, "lwr"])`  |  `r brf(P_b[, "upr"])`  | `r brf(amp_b)`  | `r brf(ca_b[1])` | `r brf(ca_b[2])`  |  
| 2      |`r brf(p1_b[, "fit"])`  | `r brf(P1_b[, "lwr"])`|  `r brf(P1_b[, "upr"])` | `r brf(amp1_b)` | `r brf(ca1_b[1])`| `r brf(ca1_b[2])` |
| 3      |`r brf(p2_b[, "fit"])`  | `r brf(P2_b[, "lwr"])`|  `r brf(P2_b[, "upr"])` | `r brf(amp2_b)` | `r brf(ca2_b[1])`| `r brf(ca2_b[2])` |
Table: Comparação dos limites do IP e CA para os modelos do Ex. 2.

### Mercado 3

```{r dadosEx3, echo = FALSE}
set.seed(seed)
area <- c(360, 360, 360, 360, 480, 480, 360, 360, 480, 480)
situacao <- c(rep(0, 4), rep(1, 2), rep(0, 4))
VU <- 5000 - 5*area + 780*situacao + rnorm(10, 0, epsilonsd)
dados <- data.frame(VU, Area = area, 
                    Situacao = factor(situacao, 
                                      levels = c(0, 1), 
                                      labels = c("Meio de Quadra", "Esquina"))
)
```

Neste terceiro caso, foram seguidos exatamente os mesmos passos dos casos
anteriores. Contudo, para a geração dos valores, foi utilizado um peso maior
para os lotes de esquina, conforme a equação \ref{eq:ex3}, onde 
$\varepsilon \sim \mathcal{N}(0, \,100^2)$ :

\begin{equation} \label{eq:ex3}
VU = 5000 - 5 \cdot area + 780 \cdot situacao + \varepsilon
\end{equation}

```{r}
xtable(dados, 
       digits = 0,
       label = "tab:ex3",
       caption = "Dados para o M. 3.",
       align = rep("c", 4))  %>%
  xtable2kable(include.rownames = FALSE)  %>% 
  kable_styling(latex_options = "striped", 
                full_width = F, 
                position = "float_right")
```

Os dados gerados encontram-se na tabela \ref{tab:ex3}. Deve-se reparar que,
aqui, erros à parte, foi adicionado R$780,00 ao valor de um lote em situação de
esquina em relação a um lote em situação de meio de quadra, o que significa um
adicional de exatamente 30% ao valor do lote-padrão de 480 $m^2$ em situação de
meio de quadra, que apresenta valor 'real' de R\$2.600,00. Ou seja, o valor do
lote que se pretende prever, de 480 $m^2$, em situação de esquina, neste
mercado, terá valor 'real' de `r reais()(5000 - 5*480 + 780)`. Novamente, então, 
foram ajustados três modelos, como os descritos no exemplo anterior. 

Na tabela \ref{tab:tab3} são mostrados os coeficientes dos três
modelos ajustados. Novamente fica claro que: no primeiro modelo (coluna 1), os
coeficientes ajustados estão viesados, no segundo e no terceiro modelos eles
estão estimados corretamente, porém o terceiro modelo é mais preciso (menores
intervalos de confiança).

```{r, echo = FALSE}
# fit_c <- lmrob(VU ~ Area, dados, x = TRUE, y = TRUE)
fit_c <- lm(VU ~ Area, dados, x = TRUE, y = TRUE)
fit1_c <- lm(VU ~ Area, dados, subset = -c(5,6), x = TRUE, y = TRUE)
fit2_c <- lm(VU ~ Area + Situacao, dados, x = TRUE, y = TRUE)
```

```{r}
# Adjust standard errors
cov1         <- vcovHC(fit_c, type = "HC3")
robust_se    <- sqrt(diag(cov1))
```

```{r, results='asis', echo = FALSE}
stargazer::stargazer(fit_c, fit1_c, fit2_c, label = "tab:tab3", 
                     title = "Comparação de modelos para o Ex. 3.",
                     type = type, header = FALSE, table.placement = "H",
                     ci = TRUE, ci.level = .80, report = "vcst*", single.row = T,
                     intercept.bottom = FALSE, intercept.top = TRUE,
                     decimal.mark = ",", digit.separator = ".",
                     digits = 1, star.cutoffs = c(0.30, 0.20, 0.10),
                     se = list(robust_se, NULL, NULL),
                     notes = c("Reportados erros robustos para modelo da coluna 1."),
                     notes.label = "Notas:", font.size = "footnotesize")
# tab_model(fit_c, fit1_c, fit2_c,
#           dv.labels = c("Modelo 1", "Modelo 2", "Modelo 3"),
#           show.ci = 0.80,
#           show.stat = TRUE,
#           show.intercept = TRUE)
```

```{r}
newx <- seq(min(dados$Area), max(dados$Area), by=2.5)
pred.int3 <- predict(fit_c, newdata = data.frame(Area = newx), 
                         interval = "prediction", level = 0.80)
pred.int3 <- data.frame(VU = pred.int3, Area = newx)
```

```{r figex3, fig.cap="Gráfico dos modelos para o terceiro mercado."}
plt3 <- ggplot(pred.int3, aes(x = Area, y = VU.fit)) +
  geom_point(data = dados, aes(x = Area, y = VU, color = Situacao), size = .5) +
  geom_smooth(data = pred.int3, aes(ymin = VU.lwr, ymax = VU.upr), 
              stat = "identity", size = .5, color = "orange") +
  geom_abline(intercept = coef(fit1_c)[1], slope = coef(fit1_c)[2], 
              color = "blue", size = .5, lwd = 1) +
  geom_segment(aes(x = 360, y = 1.15*pred.int3[1, 1], 
                   xend = 480, yend = 1.15*pred.int3[49, 1]),
               color = "orange", lty = "ff") +
  geom_segment(aes(x = 360, y = 0.85*pred.int3[1, 1], 
                   xend = 480, yend = 0.85*pred.int3[49, 1]),
               color = "orange", lty = "ff") +
  scale_y_continuous(name = "VU", expand = c(0.025, 0.025), limits = c(2000, 4500)) +
  scale_x_continuous(expand = c(0.025, 0.025), limits = c(360, 480)) +
  theme(legend.position = 'bottom', legend.direction = "horizontal",
        axis.text.x = element_text(color = "grey20", size = 8),
        axis.text.y = element_text(color = "grey20", size = 8))
```

```{r, echo = FALSE}
p_c <- Predict(fit_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
               interval = "confidence", level = 0.80, vcov. = hccm(fit_c))
p1_c <- predict(fit1_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
p2_c <- predict(fit2_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
P_c <- Predict(fit_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80, vcov. = hccm(fit_c))
P1_c <- predict(fit1_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
P2_c <- predict(fit2_c, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
amp_c <- amplitude(P_c)
amp1_c <- amplitude(P1_c)
amp2_c <- amplitude(P2_c)
ca_c <- campo_arbitrio(p_c)
ca1_c <- campo_arbitrio(p1_c)
ca2_c <- campo_arbitrio(p2_c)
```

A presença de heteroscedasticidade mais uma vez pode ser verificada pelo teste
de Breusch-Pagan (não mostrado), ou também através da análise da Figura
\ref{fig:hetero}: para oslotes de maior área (e menor VU), a grande discrepância
entre os resíduos dos lotes em situação de esquina (resíduos positivos, pontos 5
e 6) e dos resíduos dos lotes em situação de meio de quadra (resíduos negativos
de menor área) fazem com que a variância dos erros naquele ponto sejam maiores
do que a variância dos lotes de maior área, causando a heteroscedasticidade.

```{r hetero, fig.cap = "Gráfico de resíduos vs. valores ajustados."}
plot(fit_c, which = 1, id.n = 4)
```

Ou seja, a análise dos resíduos do modelo indica uma estrutura presente nestes
resíduos. Como se sabe, esta estrutura ali presente deve-se à variável omitida:
como o modelo foi ajustado considerando dados de diferentes características sem
uma variável que represente esta mudança de característica, o valor do
coeficiente da variável `Area` tornou-se uma espécie de média ponderada entre a
redução real do valor unitário devido ao aumento da área (para os dados de
meio de quadra) e a redução menor dos valores unitários que ocorre nos dados em
situação de esquina. Os erros, portanto, apesar de simétricos (há exatamente 2
dados de cada característica com 480 $m2$) serão maiores, em magnitude, para os
imóveis de 480 $m^2$ do que para os imóveis de 360 $m^2$, donde advém a
heteroscedasticidade. Na segunda abordagem, como os dados em situação de esquina
foram removidos junto com a variável, este modelo resulta homoscedástico (não 
mostrado).

Os intervalos de confiança e predição \@80% para os três modelos, bem como os
limites do CA podem ser vistos na tabela 9 [^3] abaixo:


| Modelo | Estimativa central     |$IP_{inf}$             | $IP_{sup}$              | Amplitude IP (%)    | $CA_{inf}$       | $CA_{sup}$        |
|:-------|-----------------------:|----------------------:|------------------------:|----------------:|-----------------:|------------------:|
| 1      |`r brf(p_c[, "fit"])`   |`r brf(P_c[, "lwr"])`  |  `r brf(P_c[, "upr"])`  | `r brf(amp_c)`  | `r brf(ca_c[1])` | `r brf(ca_c[2])`  |  
| 2      |`r brf(p1_c[, "fit"])`  | `r brf(P1_c[, "lwr"])`|  `r brf(P1_c[, "upr"])` | `r brf(amp1_c)` | `r brf(ca1_c[1])`| `r brf(ca1_c[2])` |
| 3      |`r brf(p2_c[, "fit"])`  | `r brf(P2_c[, "lwr"])`|  `r brf(P2_c[, "upr"])` | `r brf(amp2_c)` | `r brf(ca2_c[1])`| `r brf(ca2_c[2])` |
Table: Comparação dos limites do IP e CA para os modelos do Ex. 3.

[^3]: O intervalos de predição para o modelo da primeira abordagem foram 
calculados com a utilização de erros robustos.

Deve-se notar que, com a primeira abordagem, como o modelo é viesado, o campo de
arbítrio ainda resultou suficiente para absorver os efeitos da variável omitida,
o que não acontece com a segunda abordagem. Para o terceiro modelo, novamente,
as previsões centrais são boas e o IP é pequeno.

\newpage

### Mercado 4

```{r dadosEx4, echo = FALSE}
set.seed(seed)
area <- c(360, 360, 360, 360, 480, 480, 360, 360, 480, 480)
situacao <- c(rep(0, 4), rep(1, 2), rep(0, 4))
VU <- 5000 - 5*area + 1500*situacao + rnorm(10, 0, epsilonsd)
dados <- data.frame(VU, Area = area, 
                    Situacao = factor(situacao, 
                                      levels = c(0, 1), 
                                      labels = c("Meio de Quadra", "Esquina"))
)
```

```{r}
xtable(dados, 
       digits = 0,
       label = "tab:ex4",
       caption = "Dados para o M. 4.",
       align = rep("c", 4))  %>%
  xtable2kable(include.rownames = FALSE)  %>% 
  kable_styling(latex_options = "striped", 
                full_width = F, 
                position = "float_right")
```

Neste quarto caso, foram seguidos exatamente os mesmos passos dos casos
anteriores. Contudo, para a geração dos valores, foi adotada uma majoração maior
para os lotes em situação de esquina, conforme \ref{eq:ex4}, onde 
$\varepsilon \sim \mathcal{N}(0, \,100^2)$ . O que se buscou foi uma situação em 
que a amplitude do IP para a primeira abordagem extrapole a amplitude do CA.

\begin{equation} \label{eq:ex4}
VU = 5000 - 5 \cdot Area + 1500 \cdot Situacao + \varepsilon
\end{equation}

Os dados gerados encontram-se na tabela \ref{tab:ex4}. Deve-se notar que os
valores das estimativas centrais para os lotes em estudo, de 480 $m^2$, em
situação de meio de quadra ou esquina, com este modelo, serão de  
`r reais()(5000 - 5*480)` e `r reais()(5000 - 5*480 + 1500)`, respectivamente. O
que se pretende é simular se num mercado com uma maior variabilidade o campo
de arbítrio pode ser uma bom fator limitante. Novamente, então, foram ajustados 
três modelos, como os descritos no exemplo anterior.

Na tabela \ref{tab:tab4} são mostrados os coeficientes dos três modelos
ajustados. Novamente fica claro que: no primeiro modelo (coluna 1), os
coeficientes ajustados estão viesados, no segundo e no terceiro modelos eles
estão estimados corretamente, porém o terceiro modelo é mais preciso (menores
intervalos de confiança).

Deve-se notar ainda na tabela \ref{tab:tab4} que o modelo da primeira abordagem,
visto na coluna (1), apresentou coeficiente de ajuste baixíssimo e coeficiente
de correlação ajustado negativo. É possível mostrar que, de acordo com o teste
de Breusch-Pagan, a hipótese da homoscedasticidade não pode ser verificada. O 
teste de significância da variável independente do modelo não foi significante 
nem para enquadrar omodelo no grau I de fundamentação da NBR 14.653-02. Ainda 
assim este modelo faz boas previsões se utilizado corretamente o IP. 

Já o modelo da segunda abordagem (coluna 2) tem um coeficiente de ajuste muito
alto e o modelo tem grau III de fundamentação, segundo a NBR 14.653-02. No 
entanto, este modelo é totalmente inadequado para fazer previsões fora dos dados
da amostra, mesmo com a utilização do CA ou do IP.

```{r, echo = FALSE}
# fit_d <- lmrob(VU ~ Area, data = dados)
fit_d <- lm(VU ~ Area, data = dados)
fit1_d <- lm(VU ~ Area, dados, subset = -c(5,6), x = TRUE, y = TRUE)
fit2_d <- lm(VU ~ Area + Situacao, dados, x = TRUE, y = TRUE)
```

```{r}
# Adjust standard errors
cov1         <- vcovHC(fit_d, type = "HC3")
robust_se    <- sqrt(diag(cov1))
```

```{r}
newx <- seq(min(dados$Area), max(dados$Area), by=2.5)
pred.int4 <- predict(fit_d, newdata = data.frame(Area = newx), 
                         interval = "prediction", level = 0.80)
pred.int4 <- data.frame(VU = pred.int4, Area = newx)
```

```{r figex4, fig.cap="Gráfico dos modelos para o quarto mercado."}
plt4 <- ggplot(pred.int4, aes(x = Area, y = VU.fit)) +
  geom_point(data = dados, aes(x = Area, y = VU, color = Situacao), size = .5) +
  geom_smooth(data = pred.int4, aes(ymin = VU.lwr, ymax = VU.upr), 
              stat = "identity", size = .5, color = "orange") +
  geom_abline(intercept = coef(fit1_d)[1], slope = coef(fit1_d)[2], 
              color = "blue", size = .5, lwd = 1) +
  geom_segment(aes(x = 360, y = 1.15*pred.int4[1, 1], 
                   xend = 480, yend = 1.15*pred.int4[49, 1]
                   ),
               color = "orange", lty = "ff") +
  geom_segment(aes(x = 360, y = 0.85*pred.int4[1, 1], 
                   xend = 480, yend = 0.85*pred.int4[49, 1]
                   ),
               color = "orange", lty = "ff") +
  scale_y_continuous(name = "VU", expand = c(0.025, 0.025), limits = c(2000, 4500)) +
  scale_x_continuous(expand = c(0.025, 0.025), limits = c(360, 480)) +
  theme(legend.position = 'bottom', legend.direction = "horizontal",
        axis.text.x = element_text(color = "grey20", size = 8),
        axis.text.y = element_text(color = "grey20", size = 8))
```


```{r, echo = FALSE}
p_d <- Predict(fit_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
               interval = "confidence", level = 0.80, vcov. = hccm(fit_d))
p1_d <- predict(fit1_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
p2_d <- predict(fit2_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "confidence", level = 0.80)
P_d <- Predict(fit_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
               interval = "prediction", level = 0.80, vcov. = hccm(fit_d))
P1_d <- predict(fit1_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
P2_d <- predict(fit2_d, newdata = data.frame(Area = 480, Situacao = "Esquina"), 
              interval = "prediction", level = 0.80)
amp_d <- amplitude(P_d)
amp1_d <- amplitude(P1_d)
amp2_d <- amplitude(P2_d)
ca_d <- campo_arbitrio(p_d)
ca1_d <- campo_arbitrio(p1_d)
ca2_d <- campo_arbitrio(p2_d)
```

O intervalo de predição \@80% para os três modelos, bem como os 
limites do CA podem ser vistos na tabela 11 [^4] abaixo:


| Modelo | Estimativa central     | $IP_{inf}$             | $IP_{sup}$             | Amplitude IP (%)   | $CA_{inf}$       | $CA_{sup}$        |
|:-------|-----------------------:|-----------------------:|-----------------------:|---------------:|-----------------:|------------------:|
| 1      |`r brf(p_d[, "fit"])`   | `r brf(P_d[, "lwr"])`  | `r brf(P_d[, "upr"])`  | `r brf(amp_d)` | `r brf(ca_d[1])` | `r brf(ca_d[2])`  |  
| 2      |`r brf(p1_d[, "fit"])`  | `r brf(P1_d[, "lwr"])` | `r brf(P1_d[, "upr"])` | `r brf(amp1_d)`| `r brf(ca1_d[1])`| `r brf(ca1_d[2])` |
| 3      |`r brf(p2_d[, "fit"])`  | `r brf(P2_d[, "lwr"])` | `r brf(P2_d[, "upr"])` | `r brf(amp2_d)`| `r brf(ca2_d[1])`| `r brf(ca2_d[2])` |
Table: Comparação dos limites do IP e CA para os modelos do Ex. 4.

[^4]: Os intervalos de predição e confiança para o modelo da primeira abordagem
foram calculados com a utilização de erros robustos.


```{r, results='asis', echo = FALSE}
stargazer::stargazer(fit_d, fit1_d, fit2_d, label = "tab:tab4", 
                     title = "Comparação de modelos para o Ex. 4.",
                     type = type, header = FALSE, table.placement = "H",
                     ci = TRUE, ci.level = .80, report = "vcst*", single.row = T,
                     intercept.bottom = FALSE, intercept.top = TRUE,
                     decimal.mark = ",", digit.separator = ".",
                     digits = 1, star.cutoffs = c(0.30, 0.20, 0.10),
                     se = list(robust_se, NULL, NULL),
                     notes = c("Reportados erros robustos para modelo da coluna 1."),
                     notes.label = "Notas:", font.size = "footnotesize")
# tab_model(fit_c, fit1_c, fit2_c,
#           dv.labels = c("Modelo 1", "Modelo 2", "Modelo 3"),
#           show.ci = 0.80,
#           show.stat = TRUE,
#           show.intercept = TRUE)
```

Deve-se notar que, neste mercado, como a majoração do valor unitário de um lote 
devido a situação esquina é grande, com a primeira abordagem o CA 
resultou insuficiente para absorver os efeitos da variável omitida. A segunda 
abordagem segue com o mesmo modelo. Já  com o terceiro modelo, novamente, as 
previsões centrais são muito boas e o IP é pequeno.

Na primeira abordagem, deve-se notar que o IP do modelo contém o valor 'real' do 
lote, enquanto que o intervalo de valores dentro das semi-amplitudes de $\pm$ 
15% do CA não contém o valor 'real' do lote, já que o valor real do lote é 
$\approx$ `r pct(P_d[, "upr"]/p_d[, "fit"]-1)` superior à estimativa do valor 
central obtida com este modelo (`r brf(p_d[, "fit"])`).

### Síntese

Em suma, para quatro mercados diferentes (4 bairros diferentes, digamos), foram
testados três tipos de abordagens diferentes: 

1. A abordagem de se considerar todos os dados, porém excluindo a variável 
`situacao`, por conta da micronumerosidade;

2. A abordagem de se desconsiderar, além das variáveis, também os dados em
situação de esquina, chegando-se assim a um modelo apenas para os lotes em
situação de meio de quadra;

3. A abordagem de utilizar-se todos os dados e todas as variáveis, apesar disso 
não ir ao encontro com o estabelecido pela norma vigente, por conta da 
micronumerosidade.

Os gráficos dos modelos para os quatro mercados podem ser vistos na Figura
\ref{fig:modelos}. Nela podem-se ver os dados do modelo acompanhados das retas
de regressão para a primeira abordagem (linha cheia, em laranja) e para as
segunda abordagem (em azul). Em cinza pode-se ver as bandas do IP para o modelo
da primeira abordagem. As linhas laranjas tracejadas representam os limites do
CA para a primeira abordagem. Deve-se notar que, como da própria definição de
IP, este possui a propriedade de conter 80% dos dados previstos para o mercado,
de maneira que, se a variável `Situacao` não está presente no modelo, o IP
deverá se alargar de maneira a conter 80% dos dados. Já num modelo ajustado com
todas as variáveis relevantes, esse intervalo é menor, já que a diferença de 
valores entre dados em situação de esquina e meio de quadra seria explicada pelo 
efeito da variável adicional e não seria absorvido pelos erros do modelo, como 
na primeira abordagem.

```{r modelos, fig.cap="Gráficos dos modelos para os quatro mercados.", out.width="100%"}
prow <- plot_grid(
  plt1 + theme(legend.position="none"), 
  plt2 + theme(legend.position="none"), 
  plt3 + theme(legend.position="none"), 
  plt4 + theme(legend.position="none"), 
  hjust = -1,
  nrow = 2,
  labels = "AUTO")
legend <- get_legend(
  # create some space to the left of the legend
  plt1 + theme(legend.box.margin = margin(0, 0, 0, 0))
)
plot_grid(prow, legend, ncol = 1, rel_heights = c(1, .1))
```

Em todos os mercados (bairros) **a terceira abordagem mostrou-se superior**,
obtendo-se um valor de estimativa central mais próxima do valor "real", e um
IP mais estreito, além de um modelo com estimação adequada para os coeficientes.

```{r}
# p <- plot_ly(data = dados, x = ~Area, y = ~Situacao, z = ~VU, opacity = 0.6) %>% 
#   add_markers()
# 
# p %>% add_surface(z = ~plane, x = ~Area, y = ~Situacao, showscale = FALSE) %>% layout(showlegend = FALSE)
#coplot(VU ~ Area|Situacao, data = dados, panel=panel.car, col="red")
```


```{r comparacao}
df <- data.frame(Mercado = c(rep("A", 3), rep("B", 3), rep("C", 3), rep("D", 3)),
                 Abordagem = rep(c(1,2,3), 4),
                 `Estimativa Central` = c(p_a[, "fit"], p1_a[, "fit"], p2_a[, "fit"],
                                          p_b[, "fit"], p1_b[, "fit"], p2_b[, "fit"],
                                          p_c[, "fit"], p1_c[, "fit"], p2_c[, "fit"],
                                          p_d[, "fit"], p1_d[, "fit"], p2_d[, "fit"]),
                 `IP superior` = c(P_a[, "upr"], P1_a[, "upr"], P2_a[, "upr"],
                                   P_b[, "upr"], P1_b[, "upr"], P2_b[, "upr"],
                                   P_c[, "upr"], P1_c[, "upr"], P2_c[, "upr"],
                                   P_d[, "upr"], P1_d[, "upr"], P2_d[, "upr"]),
                 `CA superior` = c(ca_a[2], ca1_a[2], ca2_a[2],
                                   ca_b[2], ca1_b[2], ca2_b[2],
                                   ca_c[2], ca1_c[2], ca2_c[2],
                                   ca_d[2], ca1_d[2], ca2_d[2]),
                 `Valor Real` = c(rep(2850, 3), rep(2990, 3), rep(3380, 3), rep(4100, 3))
                 )
kable(df,
      caption = "IP vs. CA em vários mercados, com diferentes abordagens.",
      booktabs = TRUE,
      digits = 1,
      format.args = list(big.mark = ".", decimal.mark = ","),
      col.names = c("Mercado", "Abordagem", "Central", "IP superior", 
                    "CA superior", "Valor 'real'")) %>%
  add_header_above(c(" ", " ", "Estimativa" = 3, " "), bold = T) %>% 
  column_spec(1, bold = T) %>%
  row_spec(0, bold = T) %>% 
  row_spec(c(1:3, 7:9) -1, extra_latex_after = "\\rowcolor{gray!6}") %>% 
  collapse_rows()
```


Com a **segunda abordagem**, obtem-se sempre o mesmo modelo, com mesmos
intervalos de predição e mesmos  campos de arbitrio, haja vista que o que se
altera nos três mercados é apenas o efeito de majoração devido à situação de
esquina, cujos dados são excluídos nesta abordagem. Apesar deste modelo obter
coeficientes ajustados corretamente para as demais variáveis, ele não é bom para
se efetuar previsões para dados com características diferentes dos dados
efetivamente utilizados na amostra. A adoção do CA  ou do IP só pode ser feita 
completamente no escuro, sem que exista nada que possa embasar a magnitude da 
majoração dos valores das estimativas centrais, que são válidas apenas para os
tipos de dados permanecentes na amostra, no caso, os dados em situação de 
meio de quadra.

Finalmente, com a **primeira abordagem**, apesar dos coeficientes serem
estimados de maneira muito pobre, viesados por conta do efeito da variável
relevante omitida, o modelo obtido é razoavelmente bom para fazer previsões,
desde que se faça o uso correto do intervalo de predição. Deve-se ter em mente
que o modelo não prevê como estimativa de valor central nem o valor para o lote
em situação de meio de quadra, nem o valor do lote em situação de esquina, mas
algo entre as duas situações. Para se prever valores dos lotes de meio de quadra
deve-se minorar os valores das estimativas centrais enquanto que para os lotes
de esquina, deve-se majorá-las[^5].

[^5]: Essentially, all models are wrong, but some are useful [@Box1986, 424].

| Mercado   | Erro $\hat Y$                   | Erro $\hat Y_{sup}$             | Erro $CA_{sup}$            |
|:---------:|--------------------------------:|--------------------------------:|---------------------------:|
| A         | `r pct(p1_a[, "fit"]/2850 - 1)` | `r pct(P1_a[, "upr"]/2850 - 1)` | `r pct(ca1_a[2]/2850 - 1)` |
| B         | `r pct(p1_b[, "fit"]/2990 - 1)` | `r pct(P1_b[, "upr"]/2990 - 1)` | `r pct(ca1_b[2]/2990 - 1)` |
| C         | `r pct(p1_c[, "fit"]/3380 - 1)` | `r pct(P1_c[, "upr"]/3380 - 1)` | `r pct(ca1_c[2]/3380 - 1)` |
| D         | `r pct(p1_d[, "fit"]/4100 - 1)` | `r pct(P1_d[, "upr"]/4100 - 1)` | `r pct(ca1_d[2]/4100 - 1)` |
Table: Erros (%) obtidos com a utilização da segunda abordagem.



Com a adoção do limite superior do IP, porém, o erro cometido com a primeira
abordagem é muito pequeno, independente do mercado estudado, como ilustra a
tabela abaixo. Os erros que se cometeriam ao se utilizar o limite superior do
CA, comum na prática da engenharia de avaliações, é muito maior do que o erro
cometido com a adoção do limite superior do IP.

| Mercado   | Erro $\hat Y$                  | Erro $\hat Y_{sup}$            | Erro $CA_{sup}$           |
|:---------:|-------------------------------:|-------------------------------:|--------------------------:|
| A         | `r pct(p_a[, "fit"]/2850 - 1)` | `r pct(P_a[, "upr"]/2850 - 1)` | `r pct(ca_a[2]/2850 - 1)` |
| B         | `r pct(p_b[, "fit"]/2990 - 1)` | `r pct(P_b[, "upr"]/2990 - 1)` | `r pct(ca_b[2]/2990 - 1)` |
| C         | `r pct(p_c[, "fit"]/3380 - 1)` | `r pct(P_c[, "upr"]/3380 - 1)` | `r pct(ca_c[2]/3380 - 1)` |
| D         | `r pct(p_d[, "fit"]/4100 - 1)` | `r pct(P_d[, "upr"]/4100 - 1)` | `r pct(ca_d[2]/4100 - 1)` |
Table: Erros (%) obtidos com a utilização da primeira abordagem.

\newpage

## Estudo de Caso 2 -- Avaliação intervalar

No estudo de caso anterior trabalhou-se com o problema da omissão de variável
relevanto no modelo. No presente estudo de caso será averiguada a previsão de
valores e avaliação de intervalos admissíveis, devido ao efeito de variáveis
não-relevantes, isto é, variáveis que agregariam pouco grau de explicação ao
modelo. Como foi visto na revisão bibliográfica, não se pode acrescentar tantas
variáveis quanto se deseje a um modelo, já que isso poderia levar a um 
superajuste do modelo, que seria bom para explicar a amostra de trabalho, mas 
não o mercado em análise. Desta forma, os (pequenos) efeitos destas variáveis 
ausentes são absorvidos pelo termo de erro, o que é desejado.

Imagine-se, então, que durante a construção de um modelo sobre uma amostra de 
mercado de apartamentos em uma zona central de uma cidade qualquer, não foram
incluídas no modelo uma série de variáveis que, apesar dos seus pequenos
efeitos, acabam por impactar de alguma forma na formação de preço dos 
apartamentos (melhor posição solar, melhor ventilação, melhor vista), já que 
o número de dados da amostra é limitado. 

Porém, um avaliador experiente, ciente das limitações do seu modelo, durante a
avaliação de um apartamento nesta zona central, se deparou com o seguinte 
problema: o imóvel avaliando possuía, para todas as características não 
contempladas no modelo, valores acima da média. Por exemplo, imagine-se que o
imóvel avaliando possuía uma vista privilegiada, uma melhor posição solar e 
ainda melhor posição em relação aos ventos, de tal maneira que, ao se avaliar 
tal imóvel pelo valor médio calculado pelo modelo estaria-se notoriamente 
sub-avaliando o mesmo. No estudo de caso anterior, a utilização do Campo de 
Arbítrio do Avaliador era possível, já que *variáveis relevantes* não foram 
*contempladas no modelo *. Neste caso, porém, o problema é diverso: o imóvel
avaliando não vale mais porque houve a omissão de uma variável relevante, mas 
sim por uma soma de pequenos efeitos que, somados, tornaram-se relevantes.

```{r dadosEx5, echo = FALSE}
set.seed(seed)
area <- rnorm(200, mean = 200, sd = .25*epsilonsd)
padrao <- rep(c(1, 2, 3, 5), each = 50)
VU <- 5000 - 5*area + 250*padrao + rnorm(200, 0, epsilonsd)
dados <- data.frame(VU,
                    Area = area,
                    Padrao = factor(padrao,
                                    levels = 1:5,
                                    labels = c("Muito Baixo", "Baixo", "Médio",
                                              "Alto", "Muito Alto"))
                    )
# dados <- data.frame(VU,
#                     Area = area,
#                     Padrao = factor(padrao,
#                                     levels = 1:5)
#                     )
```

```{r}
fit <- lm(VU ~ Area + Padrao, data = dados)
s <- summary(fit)
```

Para melhor ilustrar foram simulados dados conforme a equação \ref{eq:ex5}, onde
$Padrao_j$ é uma variável qualitativa (dicotômica em grupo) com níveis $j$ 
variando de 1 a 5 (muito baixo, baixo, médio, alto e muito alto) e 
$\varepsilon \sim \mathcal{N}(0, \,100^2)$ . Os gráficos do modelo ajustado 
com os dados simulados podem ser vistos na Figura \ref{fig:modelo}.

\begin{equation} \label{eq:ex5}
VU = 5000 - 5 \cdot Area + 250 \cdot Padrao_j + \varepsilon
\end{equation}

```{r modelo, fig.cap="Modelo ajustado para apartamentos.", fig.width = 9, out.width='70%'}
appraiseR::plotModel(fit, level = 0.80, interval = "prediction")
```

```{r}
newdata <- data.frame(Area = rep(150, 4), 
                      Padrao = c("Muito Baixo", "Baixo", "Médio", "Muito Alto"))
p <- predict(fit, newdata = newdata, interval = "confidence", level = 0.80)
P <- predict(fit, newdata = newdata, interval = "prediction", level = 0.80)
ca <- campo_arbitrio(p)
```

Imagine-se agora que o imóvel avaliando conta com área de 200 $m^2$ e padrão
médio (3). O valor da estimativa central para este imóvel, com este modelo,
seria de R$ `r brf(p[3, "fit"])` [^6].Imagine-se agora que o avaliador, convicto
que o valor do imóvel encontra-se acima do valor médio previsto para o modelo
para um imóvel com aquelas características, resolveu adotar, para o
imóvel-avaliando o valor do limite superior do intervalo de predição. Para os
dados simulados, este valor seria igual a `r brf(P[3, "upr"])`. Qual seria,
então, um intervalo razoável de valores admissíveis?

[^6]: Intervalo de confiança: [`r brf(p[3, "lwr"])`, `r brf(p[3, "upr"])`]

Como limite inferior, pode-se dizer que, muito provavelmente, o imóvel não está
abaixo da média para aquele padrão [^7]. Logo, o limite inferior do IC para a
estimativa de valor central seria um bom indicador de um limite inferior para os
valores admissíveis. E quanto ao limite superior admissível? Neste caso, seria
razoável adotar o limite superior do CA. Já que a NBR não permite que se
extrapole os 15% do CA, este é o limite superior do intervalo. Caso o limite
superior do IP ultrapasse os 15%, esta arbitragem não seria aceita, mas caso o
limite superior do IP esteja dentro do CA, é razoável que o intervalo de valores
admissíveis esteja entre o valor central e o limite superior do CA. Em suma, 
adotando-se estes critérios, o intervalor de valores admissíveis seria coerente, 
apesar de altamente assimétrico : `r brf(p[3, "lwr"])` $\leq$ 
`r brf(P[3, "upr"])` $\leq$ `r brf(ca[3, "C.A.S."])` 

[^7]: Não faz qualquer sentido, se o avaliador arbitra um valor acima da média
para o bem-avaliando, depois reportar como admissível um valor abaixo da média! 
Atualmente, isto é possível, conforme Figura A.2 da NBR 14.653-02 [-@NBR1465302].

Outra opção seria adotar um intervalo admissível sempre simétrico, com
semi-amplitude igual à menor distância entre: a) o valor adotado e o valor 
central, ou; b) o valor adotado e o limite superior/inferior do CA. Parte-se do
pressuposto que o valor arbitrado é considerado o mais provável pelo avaliador,
sendo que, em direção aos dois extremos, a probabilidade de ocorrência de valor
diminui simetricamente.Neste caso, o avaliador obteria, com este modelo, o 
seguinte intervalo de valores admissíveis: `r brf(p[3, "lwr"])` $\leq$ 
`r brf(P[3, "upr"])` $\leq$ `r brf(P[3, "upr"] + (P[3, "upr"] - p[3, "lwr"]))`  

Os dois intervalos podem ser visualizados nas áreas em azul claro na Figura 
\ref{fig:dists} A (assimétrico) e B (simétrico):

```{r dists, fig.cap="Distribuição de probabilidades dos intervalos sugeridos.", fig.height = 3, out.width="100%"}
pl1 <- 
  xpnorm(c(p[3, "lwr"], ca[3, "C.A.S."]), mean = p[3, "fit"], sd = s$sigma,
         return = "plot", system = "gg", digits=3) + 
  theme(legend.position="bottom") + 
  guides(fill=guide_legend(title="Prob."))
pl2 <- 
  xpnorm(c(p[3, "lwr"], P[3, "upr"] + (P[3, "upr"] - p[3, "lwr"])), 
         mean = p[3, "fit"], sd = s$sigma,
         return = "plot", system = "gg", digits=3) + 
  theme(legend.position="bottom") +
  guides(fill=guide_legend(title=""))
plot_grid(pl1, 
          pl2, 
          labels = "AUTO")
```


Em suma, acredita-se que uma boa melhoria em termos de normatização seria obtida
se a NBR 14.653 estabelecesse que, para o caso da arbitragem de valores (fora da
estimativa central) esta seja feita dentro do IP e do CA, simultaneamente. Além
disto, que seja reportado um intervalo de valores admissíveis coerente, conforme 
proposto acima.

Resta, porém, a solução de um outro problema: como foi visto na Figura 
\ref{fig:modelo}, os dados para este exemplo ainda foram modelados com uma 
dificuldade adicional, pois não foram inclusos, propositalmente, dados de um 
padrão em especial, a saber, dados de imóveis de padrão alto (foram excluídos os
dados simulados com esta característica na confecção do modelo, visando simular 
que os dados com esta características apresentavam micronumerosidade).

A supressão destes dados ao modelo deixa o avaliador sem condições de
utilizá-lo para fins de avaliação de imóveis com estas características [^8]. Se
o avaliador inclui os dados com aquelas características, ainda que em situação
de micronumerosidade, ele poderá fazer previsões para os imóveis com aquela
característica, ainda que menos precisas que as previsões feitas para as demais
características, com número maior de dados. É o que se nota na Figura
\ref{fig:modelo1}: para a característica 4 (padrão alto), o IP é maior do que
para as outras características, ou seja, as previsões obtidas com este modelo
para imóveis com estas características serão menos confiáveis, mas ainda assim é
melhor do que não poder obter previsão alguma. Pode ser, ainda, que as previsões
assim obtidas sejam melhores do que as que se obteriam com a interpolação numa
variável de código alocado. Mais uma vez, entende-se que seria melhor deixar o
modelo falar: se o intervalo de predição encontrado for muito grande, maior do
que o CA, descarta-se a utilização do modelo. Se for menor do que o CA, por que
não utilizá-lo?

[^8]: Seria possível modelar a variável com um código alocado e interpolar, mas
não é possível fazer isto para variáveis dicotômicas em grupo.

```{r}
dados <- rbind(
  dados, 
  data.frame(VU = c(5000, 5000), Area = c(200, 200), Padrao = c("Alto", "Alto"))
  )
fit1 <- lm(VU ~ Area + Padrao, data = dados)
s1 <- summary(fit1)
```

```{r}
newdata <- data.frame(Area = rep(150, 5), Padrao = levels(dados$Padrao))
p1 <- predict(fit1, newdata = newdata, interval = "prediction", level = 0.80)
ca <- campo_arbitrio(p1)
```

```{r modelo1, fig.cap="Modelo com dados em micronumerosidade.", fig.width = 9, out.width='70%'}
p <- appraiseR::plotModel(fit1, level = 0.80, interval = "prediction")
#p$plots + theme(legend.title=element_blank())
p
```


# Conclusão e Recomendações

Primeiramente, deve-se concluir que a micronumerosidade é um critério que
deveria ser menos rigoroso na normatização. A presença de poucos dados pode
fazer com que a estimação não seja tão boa, porém a retirada de dados ou omissão
de variáveis relevantes pode ser ainda pior. A princípio, a norma deveria deixar
o modelo dizer se a quantidade de dados é suficiente ou não: no caso da variável
apresentar significância estatística, por que removê-la? Além do mais, dados 
preciosos podem estar sendo jogados fora (no exemplo, 20% deles), o que faz com 
que a estimação dos outros coeficientes do modelo piorar.

A opção de apenas remover a variável, mantendo-se os dados, pode ser feita, mas
com o devido cuidado: o(s) valor(es) do(s) coeficiente(s) estimados para a(s)
variável(eis) efetivamente utilizada(s) deve(m) ser visto(s) com ressalvas, pois
nele(s) estão embutidos os efeitos da variável omitida. Além do mais, para que
isto seja possível, é necessário rever outros itens da NBR 14.653-02,
relacionados aos desvios em relação às hipóteses da inferência clássica: com a
omissão de variável relevante, conforme demonstrado, o esperado é que algumas
hipóteses sejam violadas, mas isto não é impeditivo da utilização do modelo: a
utilização de erros robustos (Eicker-White) para estimar os intervalos de
confiança e predição podem solucionar o problema a contento [ver @zonato2018]. É
necessário, contudo, normatizar este procedimento.

Com relação ao Campo de Arbítrio do Avaliador, a princípio entende-se que este
deva ser mantido na NBR 14.653-02, até por compatilidade com as outras partes da
NBR 14.653, mas também por compatibilidade com os outros métodos, como os
métodos de avaliação por fatores. A sua aplicação, contudo, deveria ser revista:
primeiramente, entende-se que é um erro a livre utilização de todo o intervalo
do CA sem um critério que possa dar embasamento a esta utilização. A amplitude
máxima do CA deveria ser mantida, porém como limitante superior/inferior para os
intervalos de valores admissíveis, não como uma opção de livre escolha para o
avaliador arbitrar qualquer valor dentro deste intervalo, como ocorre na atual
normatização. Entende-se que seria boa regra permitir ao avaliador que arbire
valores limitados simultaneamente pelo CA e pelo IP, utilizando os limites
destes intervalos como limites do intervalo de valores admissíveis, deixando
claro que a probabilidade de ocorrência diminui (simetricamente) à medida que se
afasta do valor arbitrado.

# Referências {-}