12-Regressao-Multipla.Rmd

---
editor_options: 
  markdown: 
    wrap: 72
---

# - Regressao Multipla

```{r package loads 12, echo=TRUE, message=FALSE, warning=FALSE}
library(ggplot2)
library(knitr)
library(tidyverse)
library(here)
library(tidyr)
```

É muito raro que tenhamos apenas um preditor em um contexto de regressão. Talvez em um experimento de laboratório, e mesmo assim é difícil não ter outros preditores. De modo que tudo que vimos até agora é um caso particular raríssimo e pouco útil per se. Para generalizar as derivações e resultados para múltiplos preditores é conveniente utilizar álgebra linear. Basicamente, isso significa usar vetores e matrizes para representar a regressão. As derivações das fórmulas ficam bem fáceis, bem como as propriedades do modelo de regressão. Como não estamos pressupondo que vocês conheçam o básico de Álgebra Linear, iremos pular as derivações e apresentar apenas as intuições, esperando que o que vocês aprenderam com um preditor seja suficiente para entender com múltiplos preditores.

## Revisão de Matriz e Vetores

Regressão múiltipla é mais facilmente compreendida com uso de matrizes. Portanto, vamos fazer uma rápida revisão da álgebra de vetores e matrizes.

Se eu tenho um vetor $x = [ a_1, a_2, ..., a_n]$, digo que este é um vetor linha com $n$ elementos. É possível também ter um vetor coluna:

$$
\begin{align}
    x &= \begin{bmatrix}
           a_{1} \\
           a_{2} \\
           \vdots \\
           a_{n}
         \end{bmatrix}
\end{align}
$$
Eu posso somar dois vetores linhas ou dois vetores colunas, se tiverem o mesmo número de elementos. Por exemplo, dois vetores colunas.
$$
  \begin{align}
          \begin{bmatrix}
           a_{1} \\           
           \vdots \\
           a_{n}
          \end{bmatrix} +
          \begin{bmatrix}
           b_{1} \\
           \vdots \\
           b_{n}
         \end{bmatrix}
    &= \begin{bmatrix}
           a_1 + b_{1} \\           
           \vdots \\
           a_n + b_{n}
          \end{bmatrix}     
  \end{align}
$$
E posso fazer multiplicação de vetores (existem vários tipos, aqui me restringo ao produto interno ou produto ponto de vetores), desde que a gente multiplique um vetor linha por uma vetor coluna, mas não o contrário.

$$
  \begin{align}
          \begin{bmatrix}
           a_{1}, a_2, \cdots, a_{n}
          \end{bmatrix} \cdot
          \begin{bmatrix}
           b_{1} \\
           b_{2} \\
           \vdots \\
           b_{n}
         \end{bmatrix}
    &=  a_1 \cdot b_{1} + a_2 \cdot b_2 \cdots + a_n \cdot b_{n} 
  \end{align}
$$
A razão é que a multiplicação de vetores (e matrizes em geral) é basicamente multiplicar linha com coluna. No caso de um vetor coluna multiplicado por um linha, isso não é possível.

A adição e multiplicação de matrizes é basicamente a generalização da ágebra com vetores


## Modelo básico

O modelo básico de regressão linear múltipla pode ser especificado por:

1. Existem p preditores, $X_1$, $X_2$, ..., $X_p$. Não precisamos fazer suposições sobre a distribuição dos preditores, e podem ser correlacionados ou não.

2. Há uma única variável resposta, $Y$. Se houvesse mais de uma, teríamos um modelo de regressão multivariada.

3. $y_i = \alpha + \beta_1 \cdot x_{1i} + \beta_2 \cdot x_{2i} + ... + \beta_p \cdot x_{pi} + e_i$. Portanto, temos $p+1$ parâmetros ou coeficientes de regressão a estimar.

4. O erro $e_i$ possui esperança condicional zero e variância condicional constante no modelo homocedástico, e não correlacionado entre observações.

Se assumirmos normalidade do termo de erro, temos também:

5. O erro $e_i$ tem uma distribuição normal multivariada, com vetor de médias zero e matriz de variância e covariância cujos elementos fora da diagonal (covariância) são zero, e a diagonal principal é $\sigma^2$.

## Modelo com matrizes

Vejam que $\alpha + \beta_1 \cdot x_{1i} + \beta_2 \cdot x_{2i} + ... + \beta_p \cdot x_{pi} + e_i$ é uma soma de produtos, similar ao que eu tinha com vetores no exemplo acima. Exceto que $\alpha$ não multiplica nada. Então, vou considerar que tenhao um preditor cujo valor é uma constante e igual a $1$, e os demais preditores, de forma que o lado direito da equação de regressão pode ser reescrito como soma e multiplicação de matrizes. Para o caso de um preditor $y_i = \alpha + \beta_ix_{1i} + e_i$, a equação de regressão com matrizes, fica:

$$
  \begin{align}
  \begin{bmatrix}
    y_1 \\
     y_2 \\
     \vdots \\
     y_{n}
     \end{bmatrix}
     &= 
          \begin{bmatrix}
           1 & x_{11} \\
           1 & x_{12} \\
           \vdots \\
           1 & x_{1n}
         \end{bmatrix}
          \begin{bmatrix}
           \alpha \\
           \beta_1 \\
         \end{bmatrix} +
        \begin{bmatrix}
           e_1 \\
           e_2 \\
           \vdots \\
           e_n
         \end{bmatrix}
  \end{align}
$$
Se eu chamar o vetor coluna com os $y$ de $Y$, a matriz com a constante $1$ e $x_{1i}$ de $X$, o vetor de coeficientes de $B$ e o vetor de erros $\epsilon$, tenho então:

$$
Y = XB + \epsilon
$$
Veja que a generalização para $p$ preditores gera a mesma equação:

$$
  \begin{align}
  \begin{bmatrix}
    y_1 \\
     y_2 \\
     \vdots \\
     y_{n}
     \end{bmatrix}
     &= 
          \begin{bmatrix}
           1 & x_{11} & x_{21} \cdots & x_{p1}\\
           1 & x_{12} & x_{22} \cdots & x_{p2}\\
           \vdots \\
           1 & x_{1n} & x_{2n} \cdots & x_{pn}
         \end{bmatrix}
          \begin{bmatrix}
           \alpha \\
           \beta_1 \\
           \beta_2 \\
           \vdots \\
           \beta_p
         \end{bmatrix} +
        \begin{bmatrix}
           e_1 \\
           e_2 \\
           \vdots \\
           e_n
         \end{bmatrix}
  \end{align}
$$

$$
Y = XB + \epsilon
$$

A única diferença é o tamanho da matriz $X$ e $B$, em que $X$ é uma matriz $n \times (p+1)$, isto é, com $n$ linhas e $p+1$ colunas, e $B$ é uma matriz de tamanho $(p+1) \times 1$ e $Y$ e $\epsilon$ são $n \times 1$. E as suposições podem ser escritas como $\mathbb{E}[\epsilon|X] = 0$ e $\mathbb{Var}[\epsilon|X] = \sigma^2I$, em que $I$ é a matriz identidade, isto é, uma matriz cuja diagonal principal é $1$ e o resto é $zero$.

### Estimador de MQO

É possível mostrar que o estimador de mínimos quadrados ordinários é dados por:

$$
\hat{B} = (X'X)^{-1} \cdot X'Y
$$
Veja que $X'Y$ é um produto (com soma) entre $X$ e $Y$, ou seja, é como se fosse a covariância entre $X$ e $Y$, e $X'X$ se assemelha à variância de $X$. E está elevado a $-1$ porque não existe divisão em matriz, de forma que preciso multiplicar pela inversa.

## Interpretação dos coeficientes

Nosso coeficiente $\alpha$ é novamente o valor esperado do $Y$ na origem, isto é:

$$
\alpha = \mathbb{E}[Y|X_1=0, X_2=0, \cdots, X_p=0]$
$$
Em um modelo sem interações, o efeito de cada variávei $X_i$ é a contribuição separada para a resposta esperada (média). Portanto, $B_i$ mede a contribuição de como $\mathbb{E}[Y]$ muda à medida que $X_i$ (e apenas $X_i$) muda, para qualquer valor de $X_i$ (se a equação for linear nas variáveis) e para qualquer valor das demais variáveis (pressuposto de aditividade, sem interação, dos preditores).


Vamos retomar nosso modelo de previsão eleitoral e rodar no R, agora adicionando múltipals variáveis.

```{r dados-18 12, echo=TRUE, message=FALSE, cache=TRUE}
library(data.table)

# lista o nome do arquivo em csv
# unzip(here("dados", "votacao_secao_2018_BR.zip"), list = TRUE)


#read data1.csv into data frame
presid_18 <- fread(here("dados","votacao_secao_2018_BR.csv"), encoding = "Latin-1")

# Supondo que seu dataframe seja chamado df
df_resultados <- presid_18 %>%
  dplyr::filter(!NR_VOTAVEL %in% c(95,96)) %>%
  group_by(NR_ZONA, CD_MUNICIPIO, SG_UF, NR_VOTAVEL, NR_TURNO) %>%
  summarise(total_votos = sum(QT_VOTOS)) %>%
  pivot_wider(names_from = NR_TURNO, values_from = total_votos, values_fill = 0) %>%
  clean_names() %>%
  group_by(nr_zona, cd_municipio, sg_uf) %>%
  mutate(total_validos_1t = sum(x1),
         total_validos_2t = sum(x2)) %>%
  dplyr::filter(nr_votavel %in% c(13,17)) %>%
  group_by(nr_votavel) %>%
  mutate(percentual_1t = x1 /total_validos_1t,
         percentual_2t = x2 / total_validos_2t) %>%
  ungroup() %>%
  dplyr::select(-c(x1, x2, total_validos_1t, total_validos_2t)) %>%
  pivot_wider(names_from = nr_votavel, 
              values_from = c(percentual_1t, percentual_2t))

# remove
# rm(presid_18)

df_resultados %>%
  ggplot(aes(x=percentual_1t_17, y=percentual_2t_17)) + geom_point() + facet_wrap(~sg_uf) + geom_smooth(method="lm", se=F, linewidth = .5)

# modelo de regressão

reg1 <- lm(percentual_2t_17 ~ percentual_1t_17 + percentual_1t_13 + sg_uf, data = df_resultados)
summary(reg1)

# AC é a categoria de referência

```

A interpretação das variáveis, portanto, é a seguinte:

1. O intercepto mede o percentual médio no segundo turno quando todas as variáveis são zero. Ou seja, Haddad e Bolsonaro tiveram 0 pontos percentuais (não existe caso assim!) e a UF é o Acre, que é a categoria de referência.

2. A variável "percentual_1t_17" mede o efeito preditivo no voto do 2o turno do Bolsonaro de aumento de um ponto percentual no voto do primeiro turno, que é de 0,77 pontos percentuais. A variável "percentual_1t_13" mede similarmente o efeito preditivo de aumento de um ponto percentual do voto do Haddad no primeiro turno sobre o voto do Bolsonaro no 2o turno. Como esperado, a relação é negativa, isto é, quanto melhor o Haddad foi no primeiro turno, pior o Bolsonaro no segundo turno naquela seção eleitoral. E o efeito de cada UF é o efeito de estar naquela UF, em comparação com a categoria de referência, ACRE.

## Regressão múltipla versus Múltiplas Regressões Separadas ou Viés de variável omitida

Rodar uma regressão com duas variáveis (digamos), não é o mesmo que rodar duas regressões separadas, uma com cada variável. A razão é que os preditores em geral terão alguma correlação entre si. Para ver isso, suponha que o verdadeiro modelo é $Y = \alpha + \beta_1 \cdot X_1 + \beta_2 \cdot X_2 + \epsilon$. O que aconteceria se rodássemos uma regressão com um preditor apenas ($X_1$)?

Vamos usar as seguintes propriedades da covariância nessa derivação. 

1. Sejam $X$ e $Y$ duas v.a., Seja $A$ uma constante. Então, $\mathbb{Cov}[X,Y + A] = \mathbb{Cov}[X,Y] + \mathbb{Cov}[X,A]$. Como $A$ é constante, $\mathbb{Cov}[X,A] = 0$ e, portanto, $\mathbb{Cov}[X,Y + A] = \mathbb{Cov}[X,Y]$.

2. $\mathbb{Cov}[X,Y \cdot A] = A \cdot \mathbb{Cov}[X,Y]$.

Em um modelo com um único preditor, teremos: 
$Y = \alpha + \beta_1^* \cdot X_1 + \epsilon$

Designei o beta da equação com um preditor por $\beta_1^*$, para diferenciar do $\beta_1$ da verdadeira equação, com dois preditores.

Nós sabemos que $\beta_1^* = \frac{\mathbb{Cov}[X_1,Y]}{\mathbb{Var}[X_1]}$

Vamos substituir o $Y$ do modelo verdadeiro na equação do $\beta_1^*$.


\begin{align*}
\beta_1^* = \frac{\mathbb{Cov}[X_1,\alpha + \beta_1 \cdot X_1 + \beta_2 \cdot X_2 + e]}{\mathbb{Var}[X_1]} = \\
\frac{\mathbb{Cov}[X_1,\beta_1 \cdot X_1] + \mathbb{Cov}[X_1, \beta_2 \cdot X_2] + \mathbb{Cov}[X_1, e]}{\mathbb{Var}[X_1]} = \\

 \frac{\beta_1 \cdot \mathbb{Cov}[X_1, X_1] + \beta_2 \cdot \mathbb{Cov}[X_1, X_2] + \mathbb{ Cov}[X_1,e]}{\mathbb{Var}[X_1]} = \\
 
\frac{\beta_1 \cdot \mathbb{Var}[X_1] + \beta_2 \cdot \mathbb{Cov}[X_1, X_2] + 0}{\mathbb{Var}[X_1]} = \\

\frac{\beta_1 \cdot \mathbb{Var}[X_1]}{\mathbb{Var}[X_1]} + \frac{\beta_2 \cdot \mathbb{Cov}[X_1, X_2] + 0}{\mathbb{Var}[X_1]} = \\

\beta_1 + \frac{\beta_2 \cdot \mathbb{Cov}[X_1, X_2] + 0}{\mathbb{Var}[X_1]} = \\

\beta_1^* = \beta_1 + \frac{\beta_2 \cdot \mathbb{Cov}[X_1, X_2]}{\mathbb{Var}[X_1]}
\end{align*}


Vemos que a inclinação $\beta_1^*$ inclui a contribuição direta de $X_1$ via $\beta_1$ mais a contribuição indireta da correlação com $X_2$, via $\beta_2$.

Portanto, se eu rodar uma regressão com um preditor quando o verdadeiro modelo tem dois preditores, o coeficiente de $\beta_1^*$ será uma média entre $\beta_1$ e $\beta_2$. Por outro lado, se eu rodar a regressão com o modelo correto com os dois preditores, consigo que $\beta_1^*$ reflita só a contribuição de $\beta_1$.

Talvez você esteja se perguntando a essa altura: quem garante que o verdadeiro modelo possua só dois preditores? Isso é o que chamamos de viés de variável omitida. Se omitirmos da regressão uma variável $X_k$ correlacionada com $X_j$, $j \neq k$, então o coeficiente $\beta_j$ reflitirá também o efeito de $\beta_k$. 

Aqui não estamos falando de causalidade, apenas da contribuição para a previsão da nossa variável resposta. Naturalmente, antes da moderna abordagem de inferência causal por resultados potenciais de Rubin (ou redes Bayesianas em modelos estruturais de Pearl), as pessoas pensavam que, controlando para o máximo de variáveis possível, com sorte seria possível eliminar (ou reduzir a um mínimo) o viés de variável omitida e, portanto, estar seguro que $\beta_1$ estimaria o efeito causal.

Nós hoje sabemos que o modo mais seguro de pensar causalidade é usando uma das duas abordagens (as iniciadas por Rubin ou Pearl), e verificando (por exemplo com resultados potenciais) que a suposição de independência condicional (*CIA*, de Conditional Independence Assumption) é plausível para poder interpretar $\beta_1$ causalmente. Sem um modelo causal, a abordagem de introdução de regressores para controlar o viés de variável omitida não nos permite fazer inferência causal, exceto em casos muitos simples ou quando implicitamente temos garantida a validade da *CIA* (como em um experimento bem conduzido e com compliance), como é o caso das ciências naturais como física e química.


## Matriz chapéu

Uma forma interessante de visualizar as previsões do modelo é que podemos escrever $\hat{Y} = XB$. Substituindo a fórmula do $B$, temos que:
$$
\hat{Y} = X (X'X)^{-1}X'Y = \\
(X (X'X)^{-1}X')Y = \\
HY
$$
Essa equação mostra que as previsões são dadas pelas respostas observadas, ponderadas pela matriz chapéu (*hat*), $H$.

## Multicolinearidade

Até o momento não falamos sobre em que condições a matrix inversa $X'X^{-1}$ existe. Nós sabemos que nem todas as matrizes podem ser invertíveis. Matrizes com determinante zero são não-invertíveis. A intuição é como pensar que não é possível dividir um escalar por zero. O determinante é zero quando as colunas não são linearmente independentes. Ou seja, quando uma coluna (ou mais) é uma combinação linear de um ou mais colunas. No nosso caso, quando a correlação for $1$ (ou $-1$). Nesses casos, não é possível estimar os coeficientes da regressão e acontece quando temos multicolinearidade.

Softwares modernos, como R, irão "dropar" uma (ouas mais) variável(eis) se isso ocorrer, automaticamente, para evitar que a matriz não seja invertível. Assim, a menos que acorrelação seja perfeita, multicolinearidade não costuma ser um problema.

## Erro padrão Robusto

Na presença de heterocedasticidade ou correlação nos erros (como autocorrelação temporal ou autocorrelação espacial), precisamos corrigir o cálculo do erro padrão. Para explicar como é calculado o erro padrão robusto, vamos derivar o erro padrão novamente, agora com a notação matricial.

Lembremos que: $\hat{B} = (X'X)^{-1}(X'Y)$ e $Y = XB + e$. Logo, reescrevendo a equação de regressão, temos:
\begin{align}
\hat{B} = (X'X)^{-1}(X'[XB + e]) \\
 = (X'X)^{-1}(X'XB + X'e) \\
 = (X'X)^{-1}X'XB + (X'X)^{-1}X'e \\
 = B + (X'X)^{-1}X'e
\end{align}

E a partir dessa equação, podemos calcular a variância dos estimadores.


\begin{align}
\mathbb{Var}[\hat{B}|X] = \mathbb{Var}[B + (X'X)^{-1}X'e|X] \\
 = \mathbb{Var}[(X'X)^{-1}X'e|X] \\
 = (X'X)^{-1}X'\mathbb{Var}[e|X]X(X'X)^{-1} \\
 = (X'X)^{-1}X'\sigma^2IX(X'X)^{-1} \\
 = \sigma^2(X'X)^{-1}X'X(X'X)^{-1} \\
= \sigma^2(X'X)^{-1}
\end{align}

Para ficar mais familiar para a gente, posso multiplicar e dividir por $n$, que não altero a equação. Assim, temos:
$$
\mathbb{Var}[\hat{B}|X] = \frac{\sigma^2}{n}(n^{-1}(X'X)^{-1})
$$
Lembrem-e que a variância do $\hat{\beta|x}$ no modelo de regressão simples era dada por: $\frac{\sigma^2}{nS_x^2}$

Então, $\frac{\sigma^2}{n}$ é igual ao que tínhamos antes. À medida que $n$ cresce, esperamos que $X'X$ cresça, já que é uma soma sobre todos os dados $n$. Dividingo todas as entradas da matriz por $n$ compensa isso. Se a covariância amostral entre todos os preditores fossem iguais a zero (sem correlação), então quando calculássemos a inversa obteríamos apenas a variância amostral de $X$ $1/S_x^2$ na diagonal principal, e temos um termo que já conhecemos de regressão simples. Lá, como só tem um preditor, não tem como ter covariância com outro preditor.

Agora podemos falar de erro padrão robusto. Notem que em nossa derivação, a certa altura, tivemos:
$$
\mathbb{Var}[\hat{B}|X] = (X'X)^{-1}X'\sigma^2IX(X'X)^{-1}
$$
Se chamarmos $\sigma^2I$ de $\Omega$, reescrevo a equação como:

$$
\mathbb{Var}[\hat{B}|X] = (X'X)^{-1}X'\Omega X(X'X)^{-1}
$$
E essa equação, escrita desse formato, é chamada de equação sanduíche, pois temos
1. $X'\Omega X$ no meio
2. $(X'X)^{-1}$ nas pontas.
Ou seja, a carne $X'\Omega X$ vai no meio de duas fatias de pão $(X'X)^{-1}X'$.

Eu suponho que isso seja engraçado de alguma forma (ou talvez fosse uma forma menmômica de memorizar a equação). Mas o fato é que o termo sanduíche pegou de tal forma que nosso erro padrão robusto envolve trocar a "carne".

Veja que definimos $\Omega = \sigma^2I$. E $\sigma^2I$ pressupõe que temos homecedasticidade, já que a variância do erro é constante. Temos portanto de modificar a "carne" para calcular o erro padrão-robusto. Vamos fazer isso manualmente, primeiro calculando o erro padrão tradicional, e depois o que seria um erro padrão robusto.

```{r cap12-se-robusto toy, message=FALSE}
library(ggplot2)
library(tidyverse)
set.seed(123)

x <- c(1:8, 10, 15)
y <- c(5 + rnorm(8,sd = 1.5), 40, 65)
df <- data.frame(y=y, x=x)

df %>%
  ggplot(aes(x=x, y=y)) + geom_point()

fit <- lm(y ~x)
summary(fit)
```

Nós construímos um "toy model" em que não existe relação entre $x$ e $y$. Porém, por causa de dois "outliers" o R achou que existia associação entre as variáveis. Como o R calcula o erro padrão? Nossa fórmula requer $\sigma^2$, $I$ e a matriz de preditores $X$, que depois vou transpor, calcular inversa etc. No R, podemos computar cada um desses itens manualmente da seguinte forma.

```{r se-regula1, message=FALSE}
n <- length(y)
sigma2 <- sigma(fit)^2
mat_I <- diag(n)
X <- model.matrix(fit)
```

```{r se-regula2, message=FALSE}
omega <- sigma2*mat_I
bread <- solve(t(X)%*%X)
meat <- (t(X) %*% omega %*% X)
vce <- bread %*% meat %*% bread
sqrt(diag(vce))
```
 
Se inspecionarmos omega mais detalhadamente, temos:
```{r print omega}
print(omega)
```

Como falamos, é uma matriz de variância constante (homocedástica).
Em nosso "toy model", contudo, não faz sentido assumir que a variância é constante. É provável que os dois "outliers" tenham vindo de uma distribuição com variância bem maior. E isso levaria a um erro padrão maior para a o coeficiente da variável $x$. 

A questão então é como calcular essas variância diferente para esse pontos. Há inúmeras maneiras de fazer isso. A default do Stata, chamada de "HC1", em que "HC" quer dizer "Heteroscedasticity-Consistent". A fórmular para o erro padrão robusto "HC1" é a seguinte:
$$
\frac{n}{n-k}\hat{e}^2
$$
em que $\hat{e}^2$ é o quadrado dos resíduos e $k$ o número de parâmetros. Em nosso modelo simples, $k=2$. E substituirmos a variância constante por essa fórmula em nossa "carne", temos:

```{r se-robusto hc1, message=FALSE}
sigma2_hc1 <- residuals(fit)^2*(n/(n-2))
omega <- sigma2_hc1*mat_I
bread <- solve(t(X)%*%X)
meat <- (t(X) %*% omega %*% X)
vce <- bread %*% meat %*% bread
sqrt(diag(vce))
```

Vemos que o erro-padrão robusto é maior que o anterior, mas a associação ainda é significativa. O erro padrão robusto do pacote do R "sandwich" usa o "HC3", que possui outra fórmula, e tem desempenho melhor em amostras pequenas. Sua fórmula é:

$$
\frac{\hat{e}^2}{(1-h_i)^2}
$$
Aqui, $h_i$ são os valores chapéus da matriz chapéu, $H$, que vimos antes. Os valores chpéus variam entre $0$ e $1$ e quanto maior o número, mais influente a observação. Então, sabemos que nossos dois últimos números terão valores chapéu maiores, ou seja, vão inflar a variância para esses dois valores. No R, podemos calcular os valores chapéus com a função "hatvalues".

```{r se-robusto hc3, message=FALSE}
sigma2_hc3 <- residuals(fit)^2/(1 - hatvalues(fit))^2
omega <- sigma2_hc3*mat_I
bread <- solve(t(X)%*%X)
meat <- (t(X) %*% omega %*% X)
vce <- bread %*% meat %*% bread
sqrt(diag(vce))
```

E posso calcular o intervalo de confiança com o novo erro padrão. Obviamente, não iremos fazer na "mão" esse cálculo todo. Vamos usar os pacotes "sandwich" e "lmtest" para fazer o teste de hipótese com erro padrão robusto.

```{r se-sandwich hc3, message=FALSE}
library(lmtest)
library(sandwich)
coeftest(fit, vcovHC(fit, "HC3"))
```

Vemos, portanto, que nosso erro padrão é robusto a observações muito influentes que são "outliers". Uma forma de verificar se temos esse problema é ver se temos resíduos grandes e observações muito influentes (valores chapéu altos). Um comando básico do R permite inspecionar visualmente se é o caso. Pontos nos cantos direitos superiores ou inferiores indica observações exibindo influência em nosso modelo.

Resíduos grandes (ou variância não constante) pode decorrer de um modelo mal-especificado (preditor ausente, ausência de interação, efeitos não lineares etc.).

```{r high-leverage, message=FALSE}
plot(fit, which = 5)
```

Como mencionado, há vários tipos de erros padrão robustos, "HC0", "HC1", "HC2", "HC3", entre outros. Vamos falar rapidamente sobre cada uma deles. Antes disso, vamos comparar para nossa regressão como o erro padrão muda para cada um desses tipos.

```{r robust-se, , results='asis'}
library(sandwich)
library(stargazer)
library(lmtest)
m2 <- coeftest(reg1, vcovHC(reg1, type = "HC0"))
m3 <- coeftest(reg1, vcovHC(reg1, type = "HC1"))
m4 <- coeftest(reg1, vcovHC(reg1, type = "HC2"))
m5 <- coeftest(reg1, vcovHC(reg1, type = "HC3"))
stargazer(reg1, 
          m2,
          m3,
          m4,
          m5, type = "html")

# conf int
library(knitr)
coefci(reg1, vcov. = vcovHC(reg1, type = 'HC1')) %>%
  kable()
```

### Entendendo o erro padrão robusto

Vamos fazer uma simulação no R para entender os vários tipos de erro padrão robusto. Vamos criar um modelo com heterocedasticidade.

```{r robust-se-sim, results='asis'}
x <- rnorm(10)
e <- rnorm(10, 0, x^2) # o DP do erro é igual a .5*x^2
a <- 2
b <- -2
y <- a + b*x + e

df <- data.frame(y=y, x=x)

df %>%
  ggplot(aes(y=y, x=x)) + geom_point() + geom_smooth(method="lm")

m_sim <- lm(y ~ x, data=df)
# plots de checagem do modelo
library(easystats)
check_model(m_sim)

f2 <- coeftest(m_sim, vcovHC(m_sim, type = "HC0"))
f3 <- coeftest(m_sim, vcovHC(m_sim, type = "HC1"))
f4 <- coeftest(m_sim, vcovHC(m_sim, type = "HC2"))
f5 <- coeftest(m_sim, vcovHC(m_sim, type = "HC3"))
stargazer(m_sim, 
          f2,
          f3,
          f4,
          f5, type = "html")
```

```{r dados-22 12, echo=TRUE, message=FALSE, cache=TRUE}
# dados de 2022

#read  into data frame
presid_22 <- fread(here("dados","votacao_secao_2022_BR.csv"), encoding = "Latin-1")

df_resultados_22_aux <- presid_22 %>%
  dplyr::filter(!NR_VOTAVEL %in% c(95,96)) %>%
  group_by(NR_ZONA, CD_MUNICIPIO, SG_UF, NR_VOTAVEL, NR_TURNO) %>%
  summarise(total_votos = sum(QT_VOTOS)) %>%
  pivot_wider(names_from = NR_TURNO, values_from = total_votos, values_fill = 0) %>%
  clean_names() %>%
  group_by(nr_zona, cd_municipio, sg_uf) %>%
  mutate(total_validos_1t = sum(x1),
         total_validos_2t = sum(x2)) %>%
  dplyr::filter(nr_votavel %in% c(13,22)) %>%
  group_by(nr_votavel) %>%
  mutate(percentual_1t = x1/total_validos_1t,
         percentual_2t = x2/total_validos_2t) %>%
  ungroup()

df_resultados_22 <- df_resultados_22_aux%>%
  dplyr::select(-c(x1, x2, total_validos_1t, total_validos_2t)) %>%
  pivot_wider(names_from = nr_votavel, 
              values_from = c(percentual_1t, percentual_2t)) %>%
  rename(percentual_1t_17 = percentual_1t_22,
         percentual_2t_17 = percentual_2t_22)

```

Agora que importamos os dados de 22, podemos fazer nossa previsão, usando os resultados do primeiro turno.

```{r previsao 12, echo=TRUE, message=FALSE, cache=TRUE}
df_resultados_22_aux1 <- df_resultados_22_aux %>%
  ungroup() %>%
  group_by(nr_zona, cd_municipio, sg_uf) %>%
  summarise(total_validos_1t = sum(total_validos_1t))
  
#Previsão
  
previsoes <- predict(reg1, newdata = df_resultados_22, interval = "prediction", level = .95) %>%
  as.data.frame()

df_resultados_22_final <- df_resultados_22 %>%
  ungroup() %>%
  mutate(prev_perc = previsoes$fit,
         prev_perc_lower = previsoes$lwr,
         prev_perc_upper = previsoes$upr,
         validos_prev = df_resultados_22_aux1$total_validos_1t*prev_perc,
         validos_prev_lower = df_resultados_22_aux1$total_validos_1t*prev_perc_lower,
         validos_prev_upper = df_resultados_22_aux1$total_validos_1t*prev_perc_upper)

tot_valido <- sum(df_resultados_22_aux1$total_validos_1t)
df_resultados_22_final %>%
  summarise(perc_previsto = sum(validos_prev,na.rm = T)/tot_valido,
            perc_previsto_lower = sum(validos_prev_lower,na.rm = T)/tot_valido,
            perc_previsto_upper = sum(validos_prev_upper,na.rm = T)/tot_valido)

```

## Termo de interação

No trabalho de Corrêa (2015), o autor estima um modelo para explicar o efeito do desempenho econômico (crescimento do PIB e inflação) sobre o desempenho eleitoral de candidatos governistas na América Latina. E a certa altura argumenta que o efeito do crescimento do PIB e da inflação dependem (ou variam) conforme o candidato é de esquerda, centro ou de direita. Em particular, candidatos de centro seriam mais afetados por essas variáveis. Vamos simplificar aqui e considerar apenas uma variável econômica (crescimento) e tratar a variável ideologia do presidente como binária (centro ou não).

Nesse caso, dizemos que o modelo contém um termo interativo, do seguinte modo:

$$
y_i = \alpha + \beta_1 \cdot x_1 + \beta_2 \cdot x_2 + \beta_3 \cdot x_1 \cdot x_2 + e_i
$$
Nesse tipo de modelo, o efeito de $x_1$ depende do valor de $x_2$ (e vice-versa). Isso significa que temos que mudar nossa interpretação do significado dos parâmetros.

Agora, $\beta_1$ mede o efeito no $y$ médio de aumentar uma unidade $x_1$ quando $x_2 = 0$. Para ver isso, podemos substiruir $x_2 = 0$ na equação acima, e temos:

$$
y_i = \alpha + \beta_1 \cdot x_1 + \beta_2 \cdot 0 + \beta_3 \cdot x_1 \cdot 0 + e_i = \alpha + \beta_1 \cdot x_1 + e_i
$$
E nesse caso, 

\begin{align}
\mathbb{E}[Y|x_1=1] - \mathbb{E}[Y|X=0] = \\
\mathbb{E}[\alpha|x_1=1] + \mathbb{E}[\beta_1 \cdot x_1|x_1=1] + \mathbb{E}[e_i|x_1=1] - (\mathbb{E}[\alpha|x_1=0] + \mathbb{E}[\beta_1 \cdot x_1|x_1=0] + \mathbb{E}[e_i|x_1=0]) \\
= \alpha + \beta_1 - \alpha = \\
\beta_1
\end{align}

Que é o que já estamos acostumados. Porém, considere o caso em que $x_2 \neq 0$. $y_i = \alpha + \beta_1 \cdot x_1 + \beta_2 \cdot x_2 + \beta_3 \cdot x_1 \cdot x_2 + e_i$

Estamos interessados em calcular: $\mathbb{E}[Y|x_1=1] - \mathbb{E}[Y|X=0]$. Comecemos pelo primeiro termo.

\begin{align}
\mathbb{E}[Y|x_1=1] = \\
\mathbb{E}[\alpha|x_1=1] + \mathbb{E}[\beta_1 \cdot 1|x_1=1] +
\mathbb{E}[\beta_2 \cdot x_2|x_1=1] + \mathbb{E}[\beta_3 \cdot 1 \cdot x_2|x_1=1] + \mathbb{E}[e_i|x_1=1] \\
= \alpha + 1\beta_1  + \beta_2\mathbb{E}[x_2|x_1=1] + \beta_3\mathbb{E}[x_2|x_1=1] + 0 \\
= \alpha + \beta_1 + (\beta_2+\beta_3)\mathbb{E}[x_2|x_1=1]
\end{align}

Já o segundo termo, fica:
\begin{align}
\mathbb{E}[Y|x_1=0] = \\
\mathbb{E}[\alpha|x_1=0] + \mathbb{E}[\beta_1 \cdot 0|x_1=0] +
\mathbb{E}[\beta_2 \cdot x_2|x_1=0] + \mathbb{E}[\beta_3 \cdot 0 \cdot x_2|x_1=0] + \mathbb{E}[e_i|x_1=0] \\
= \alpha + 0\beta_1  + \beta_2\mathbb{E}[x_2|x_1=0] + 0\mathbb{E}[x_2|x_1=1] + 0 \\
= \alpha + \beta_2 \mathbb{E}[x_2|x_1=0]
\end{align}

Computando a diferença, temos:
\begin{align}
\mathbb{E}[Y|x_1=1] - \mathbb{E}[Y|x_1=0] = \\
\alpha + \beta_1 + (\beta_2+\beta_3)\mathbb{E}[x_2|x_1=1] -\alpha - \beta_2\mathbb{E}[x_2|x_1=0] \\
= \beta_1 + (\beta_2+\beta_3)\mathbb{E}[x_2|x_1=1] - \beta_2\mathbb{E}[x_2|x_1=0]
\end{align}

Se $x_1$ e $x_2$ forem independentes na média, então: $\mathbb{E}[x_2|x_1=1] = \mathbb{E}[x_2|x_1=0]$, de forma que a equação acima simplifica para:

\begin{align}
\mathbb{E}[Y|x_1=1] - \mathbb{E}[Y|x_1=0] = \beta_1 + \beta_3\mathbb{E}[x_2|x_1=1]
\end{align}

Portanto, a interpretação do que é o efeito marignal ou preditivo Em um modelo com termo de interação é mais complexa do que em uma regressão sem termo de interação.

### Variância do termo de interação

O erro-padrão do $\beta$ é dado por: $\mathbb{Var}[\beta|X] = \sigma^2 (X'X)^{-1}$, que podemos reescrever para: $\sigma^2(X'X)^{-1} = \frac{\sigma^2}{n}n^{-1}(X'X)^{-1}$.
Lebrem-se que para um único preditor, $\mathbb{Var}[\beta|X] = \frac{\sigma^2}{nS^2_x}$

Vejam que $\frac{\sigma^2}{n}$ é o que tínhamos antes. E $n^{-1}(X'X)^{-1}$ é uma matriz de variância e covariância. Se as covariâncias forem zero, na diagonal principal temos as variâncias amostrais de cada preditor, de tal modo que:

Supondo que os preditores \(X_1, X_2, \ldots, X_p\) são não correlacionados, a matriz \(X'X\) será aproximadamente diagonal. Assim, podemos representar a matriz de variância e covariância apenas com as diagonais principais, onde cada elemento \(S_{x_j}^2\) na diagonal representa a variância amostral do preditor \(X_j\). Dessa forma, temos:

\[
\mathbb{Var}[\hat{\beta} \mid X] = \frac{\sigma^2}{n} 
\begin{bmatrix}
\frac{1}{S_{x_1}^2} & 0 & \cdots & 0 \\
0 & \frac{1}{S_{x_2}^2} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \frac{1}{S_{x_p}^2} \\
\end{bmatrix}.
\]

No caso de um único preditor (\(p = 1\)), essa expressão se reduz para a fórmula clássica:
\[
\mathbb{Var}[\hat{\beta} \mid X] = \frac{\sigma^2}{n S_x^2},
\]
onde \(S_x^2\) é a variância amostral do preditor \(X\). Assim, cada elemento da diagonal principal da matriz \(\mathbb{Var}[\hat{\beta} \mid X]\) representa a variância de \(\hat{\beta}_j\) para o preditor \(X_j\), dada por \(\frac{\sigma^2}{n S_{x_j}^2}\), generalizando a fórmula da variância para \(p\) preditores.

Obviamente, se houver correlação entre preditores, a fórmula se torna mais complicada do que a original. No caso de termo de interação, o pressuposto é de que há correlação, então fiquemos com a fórmula matricial, que é mais simples:

$$
\mathbb{Var}[\beta|X] = \sigma^2 (X'X)^{-1}
$$
E em geral a variância de cada coeficiente depende dos demais. E quanto maior a covariância entre duqas variáveis, maior a variância (e portanto, maior o erro-padrão).

O que importa para nós aqui é na verdade o seguinte.
Sejam $X$ e $Y$ duas variáveis aleatórias dependentes (cuja covariância não é zero). Então, $Var(X+Y) = Var(X) + Var(Y) + 2*Cov(X,Y)$. Portanto, no caso de nosso modelo com dois regressores e termo de interação, temos:

$$
Var(\beta_1 + \beta3 x_2|X) = Var(\beta_1|X) + Var(\beta3 x_2|X) + 2Cov(\beta_1, \beta3) = Var(\beta_1|X) + x_2^2Var(\beta3|X) + 2x_2Cov(\beta_1, \beta3|X)
$$
Como as variâncias são os erros-padrão ao quadrado, podemos computar manualmente a variância do efeito marginal, e depois calcular o erro-padrão. No R, porém, usualmente usamos alguns dos pacotes que permite a visualização desses efeitos já computados corretamente.


```{r interaction-term}
library(marginaleffects)
library(sjPlot)
set.seed(123)
n <- 100
x1 <- rnorm(n)
x2 <- rbinom(n, size=1, p=.5)
a <- 2
b1 <- 1.5
b2 <- -1
b3 <- 1
y <- a+ b1*x1 + b2*x2 + b3*x2*x1 + rnorm(n)

reg1 <- lm(y ~ x1*x2)
summary(reg1)

# Pacote marginaleffects
plot_slopes(reg1, variables = "x1", condition = "x2")

# pacote sjPlot
plot_model(reg1, type = "pred", terms = c("x1", "x2"))

# Truque do Mcdermont
reg2 <- lm(y ~ factor(x2) / x1) # lm(y ~ x2 + x2:x1)
summary(reg2)
```


 
### Calculando o erro padrão do termo de interação



## Análise de sensibilidade causal

Diante do problema de viés de variável omitida, podemos realizar análises de sensibilidade. Essa é uma área que em que houve bastante desenvolvimentos recentes. Apresento aqui o trabalho do Cinelli e Hazlet (2019). Antes, porém, vamos revisar o R-quadrado e derivar algumas formas de definir essa quantidade.

### Derivação do R-quadrado

Será útil reescrever o viés de variável omitida em termos do R-quadrado. Então, vamos derivar o R-quadrado.
Antes vamos definir algumas quantidades que serão úteis. Seja $T$ quanto o vetor de respostas $Y$ se desvia da média amostral $\bar{Y}$.
$$
T = Y - \bar{Y}
$$
E $T$ pode ser decomposto em duas partes. Primeiro, a porção do desvio em $T$ que é determinado pelo modelo de regressão, dada pelo vetor $\hat{M}$, definido como:
$$
\hat{M} = \hat{Y} - \bar{Y}
$$

E em segundo lugar, a porção do desvio em $T$ determinado pelos resíduos:
$$
\hat{e} = Y - \hat{Y}
$$
Portanto, $T = \hat{M} + \hat{e}$.
Se nós somarmos o quadrado dessas três quantidades para toda a amostra, teremos:
Soma dos Quadrados Totais (em inglês SST)
$$
SST = \sum{(y_i - \bar{y_i})^2}
$$

Soma dos Quadrados Médios (SQM)
$$
SSM = \sum{(\hat{y_i} - \bar{y})^2}
$$

Soma dos Quadrados dos Resíduos (SSE)
$$
SSE = \sum{\hat{e_i}^2} 
$$

Vamos mostrar que, do mesmo jeito que $\hat{T} = \hat{M} + \hat{E}$, $SST = SSM + SSE$

Lembremos que $\sum{\hat{e_i}} = 0$ e $\sum{\hat{e_i}\hat{y_i}} = 0$
\begin{align}
SST = \sum{(y_i - \bar{y_i})^2} = \\
\sum{(\hat{e_i} + \hat{y_i} - \bar{y_i})^2} = \\
\sum{(\hat{e_i} + (\hat{y_i} - \bar{y_i}))^2} = \\
\sum{\hat{e_i}^2 + 2\hat{e_i}(\hat{y_i} - \bar{y_i}) + (\hat{y_i} - \bar{y_i})^2} = \\
\sum{\hat{e_i}^2} + 2\sum{\hat{e_i}(\hat{y_i} - \bar{y_i})} + \sum{(\hat{y_i} - \bar{y_i})^2} = \\
SSE + 2\sum{(\hat{e_i}\hat{y_i} - \hat{e_i}\bar{y_i})} + SSM = \\
SSE + 2\sum{(\hat{e_i}\hat{y_i}} - 2\sum{\hat{e_i}\bar{y_i})} + SSM = \\
SSE + 0 - 0 + SSM = \\
SST = SSE + SSM = \\
\end{align}

A intuição dessas quantidades é que SST é uma medida da variação na resposta, $Y$ (é a variância sem dividir por $n$); SSM é uma medida dessa variação em $Y$ preditiva ou captada pela variação do modelo, e SSE a variação em $Y$ não capturada pelo modelo. Note também que se eu dividir ambos os lados da equação por $n$, tenho a variância do lado esquerdo, e não altero a igualdade da equação.

E o $R^2$ é definido pela razão entre SSM e SST:
$$
R^2 = \frac{SSM}{SST} = \frac{Var(\hat{y})}{Var(y)}
$$
O $R^2$ pode ser reescrito de várias maneiras. 


Veja que $y = \hat{y} + \hat{e}$, logo, $\hat{y} = y - \hat{e}$. E note que $Cov(\hat{y}, \hat{e}) = 0$. Portanto, temos:

$$
R^2 = \frac{Var(y - \hat{e})}{Var(y)} = 1 - \frac{Var(\hat{e})}{Var(y)} = 1 - \frac{Var(Y^{\perp X}))}{Var(y)}
$$
Ou seja, o R-quadrado é $1$ menos a variância em $Y$ não é predita por $X$, dividido pela variância de $Y$. Vamos mostrar agora que o R-quadrado pode ser também definido como a correlação (ao quadrado) entre $Y$ e $\hat{y}$.

Primeiro, vamos mostrar que
$$
R^2 = \frac{Cov(Y, \hat{y})}{Var(y)}
$$

$$
Cov(Y, \hat{y}) = Cov(\hat{y} + \hat{e}, \hat{y}) = Var(\hat{y}) + Cov(\hat{e}, \hat{y}) = Var(\hat{y}) 
$$
E podemos mostrar que:
$$
Var(\hat{y}) = Var(\hat{\alpha} + \hat{\beta} X) = Var(\hat{\beta}X) = \hat{\beta}^2 Var( X)
$$
Portanto,
$$
R^2 = \frac{\hat{\beta}^2 Var(X)}{Var(y)}
$$
E a fórmula do $\hat{\beta} = Cov(X,Y)/Var(x)$. Logo:

$$
R^2 = \frac{ Cov(X,Y)^2 Var(X)}{Var(x)^2 Var(y)} = \frac{ Cov(X,Y)^2}{Var(x)Var(y)}
$$
Designando a variância de $x$ por $S_x^2$ e similarmente para $y$, temos:

$$
R^2 = \frac{ Cov(X,Y)^2}{S_x^2S_y^2} = \left(\frac{ Cov(X,Y)}{S_xS_y}\right)^2
$$
Vamos agora definir o R-quadrado parcial, para regressão múltiplas.

Seja uma regressão completa (full) $y_i = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + e_i$ e a regressão restrita, dada por: $y_i = \alpha + \beta_1 x_{1i}  + e_i$. O R-quadrado parcial é definido como:
$$
R^2_{partial} = \frac{SSE(reduced) - SSE(full)}{SSE(reduced)}
$$
Veja que posso definir o R-quadrado parcial a partir da regressão restria só com intercepto, sem preditor, e a regressão completa com um preditor. No modelo restrito, $\hat{y} = \bar{y}$, o que implica que $SSE(reduced) = SST$. Logo, temos:

$$
R^2_{partial} = \frac{SST - SSE(full)}{SST} = 1 - \frac{SSE}{SST}
$$
Que é nossa definição de R-quadrado em modelo de regressão simples.

Então, o R-quadrado parcial me dá uma medida da proporção da variância predita pelo modelo completo que não pode ser predita pelo modelo simples.Ainda em outras palavras, o R quadrado parcial mede o poder preditivo adicional do modelo reduzido comparado ao completo. No caso de modelos com uma e duas variáveis, estou vendo o poder preditivo de adicionar uma variável ao modelo. Em termos de notação, podemos escrever o R-quadrado parcial do seguinte modo:

$$
R^2_{y \sim x1|x2} = 1 - \frac{Var(y^{\perp x_1, x_2})}{Var(y^{\perp x_1})} = \left(Cor(x_1^{\perp x_2}, y^{\perp x_2})\right)^2
$$

```{r cap12-r-quadrado-parcial, echo=TRUE, message=FALSE}
library(here)
library(data.table)
library(tidyverse)
library(sjlabelled) # pra remover labelled variables
library(haven)
library(janitor)
library(lubridate)

## dados
# https://www.latinobarometro.org/latContents.jsp

lat_bar23 <- sjlabelled::read_spss(here("Dados", "Latinobarometro_2023_Eng_Spss_v1_0.sav"),
                               drop.labels = TRUE) 

lat_bar23 <- lat_bar23 %>%
  mutate(S17 = as.Date(as.character(S17), format="%Y%m%d")) 

lat_bar23 <- lat_bar23 %>%
  janitor::clean_names()

# get_label(lat_bar23)

lat_bar23 <- lat_bar23 %>%
  mutate(data_base = as.Date(paste(diareal, mesreal, "2023", sep="-"), "%d-%m-%Y"),
         idade = year(as.period(interval(s17,data_base))),
         econ_12_meses = ifelse(p6stgbs %in% c(1,2), "better", 
                                ifelse(p6stgbs == 8, NA, "other")),
         econ_12_meses = relevel(as.factor(econ_12_meses), ref = "other"),
         aprovacao_presidente = ifelse(p15stgbs == 0, NA, p15stgbs),
         ideologia = ifelse(p16st %in% c(97, 98, 99), NA, p16st),
         votaria_governo = ifelse(perpart == 4, NA,
                                  ifelse(perpart == 1, 1, 0)),
         genero = factor(sexo, labels = c("homem", "mulher")),
         evangelico = ifelse(s1 %in% c(0,98), NA,
                             ifelse(s1 %in% c(2,3,4,5), 1, 0))) # não considera adventista, testemunha Jeová, Mórmon
br_latbar_23 <- lat_bar23 %>%
  mutate(idenpa = remove_all_labels(idenpa)) %>% # haven_labelled problems
  filter(idenpa == 76) %>% ## seelciona brasil
  filter(!is.na(votaria_governo) & !is.na(evangelico) & !is.na(ideologia) & !is.na(econ_12_meses))

glimpse(br_latbar_23)

reg_full <- lm(votaria_governo ~ ideologia + idade + genero + econ_12_meses + evangelico, data=br_latbar_23 )
summary(reg_full)
residuos_full <- resid(reg_full)

reg_res_evan <- lm(votaria_governo ~ ideologia + idade + genero + econ_12_meses, data=br_latbar_23 )
summary(reg_res_evan)
residuos_res_evan <- resid(reg_res_evan)

partial_R2_evangelico <- (sum(residuos_res_evan^2) - sum(residuos_full^2))/sum(residuos_res_evan^2)
print(partial_R2_evangelico)

# alternative
partial_R2_evangelico_alt <-  1 - var(resid(reg_full))/var(resid(reg_res_evan))
print(partial_R2_evangelico_alt)

# another alternative
reg_evan <- lm(evangelico ~ ideologia + idade + genero + econ_12_meses , data=br_latbar_23)

partial_R2_evangelico_alt1 <- cor(resid(reg_evan), resid(reg_res_evan))^2
print(partial_R2_evangelico_alt1)

 ## Cohen partial f2
r_squared_full <- summary(reg_full)$r.squared
r_squared_res <- summary(reg_res_evan)$r.squared

partial_f2 <- (r_squared_full - r_squared_res)/(1-r_squared_full)
print(partial_f2)

```
E agora podemos reescrever o viés de variável omitida em termos do R-quadrado parcial.

### Viés de variável omitida

Nós vimos que o viés de variável omitida, em uma amostra, pode ser explicitado por:
$$
\hat{\beta_1^*} = \hat{\beta_1} + \frac{\hat{\beta_2} \cdot \mathbb{Cov}[X_1, X_2]}{\mathbb{Var}[X_1]}
$$
Vamos reescrever essa equação do seguinte modo. Seja $\frac{\mathbb{Cov}[X_1, X_2]}{\mathbb{Var}[X_1]} = \hat{\gamma}$, então noss equação fica:
$$
\hat{\beta_1^*} = \hat{\beta_1} + \hat{\beta_2} \cdot \hat{\gamma}
$$

Em que $\hat{\beta_2}$ mede o "impacto" da variável omitida sobre Y. Ou seja, ela mede a diferença na esperança linear da resposta (diferença média) entre indivíduos que diferem em uma unidade na variável omitida e têm o mesmo status no tratamento, bem como o mesmo valor para todas as variáveis de controle restante. Digamos que estou interessado em estimar o efeito causal da experiência (variávle binária, experiente ou não experiente) sobre votação. Candidatas mais experientes devem ter mais votos, mas também maior visibilidade e reconhecimento, o que viesa a regressão (efeito da experiência). Suponha que controlei para arrecadação de campanha. Então, $\hat{\gamma}$ mede a diferença esperada (média) na votação entre pessoas experinetes (tratamento = 1) com uma mesma dada arrecadação de campanha, mas que difiram em uma unidade em visibilidade. 

Já a Quantidade $\hat{\gamma}$ parece que mede o impacto de Z sobre o tratamento, mas na verdade ela é resultado da regressão reversa, isto é, $Z = \hat{\alpha_1} + \hat{\gamma} T + \hat{\psi} X + \hat{u}$. Ou seja, mede a diferença média na variável de confusão entre indivíduos que diferem de uma unidade no tratamento. Em nosso exemplo, como candidatas experientes diferem em visibilidade, em média, de não experientes. Ainda em outras palavras, mede o "balanceamento" (condicional aos controles) entre indivíduos alocados pro tratamento e controle no que diz respeito à variável omitida.

Essa formulação é útil para pensar sobre como o viés de variável omitida impacta nossas estimativas. Quanto maior o "imbalance", maior será $\hat{\gamma}$. Similarmente, quanto maior o impacto da variável omitida na resposta, maior $\hat{\beta_2}$. Por isso que é comum checar o imbalance nas variáveis de controle. Se ele é alto, pode significar alto viés das variáveis omitidas, se o imbalance delas for correlacionado com o imbalance das variáveis de controle.

Numa regressão mais completa, da forma, $y = \alpha + \delta T + \beta X + \beta_2 Z + e$, em que $T$ é o tratamento, $X$ controles e $Z$ a variável omitida, o viés pode ser reescrito como:
$$
\hat{<viés>} = \hat{\beta_2} \cdot \hat{\gamma}
$$
Em que $ hat{\gamma} = \frac{Cov(T^{\perp X},Z^{\perp X} )}{Var(T^{\perp X})}$

### Viés de variável omitida em função do R-quadrado

$$
\hat{<viés>} = \sqrt{\frac{R^2_{y\sim Z|T,X} R^2_{T \sim Z|X}}{1 - R^2_{T \sim Z|X}}} \cdot \frac{sd(Y^{\perp X,T})}{sd(T^{\perp X})}
$$

Vamos olhar cada componente dessa parametrização do viés:

1. $R^2_{y\sim Z|T,X}$ é o R-quadrado parcial de $Z$ com relação à variável resposta $y$ controlando para $T$, o tratamento, e $X$ os controles. Ou seja, quanto $Z$ contribui para explicar a variação em $Y$, após descontar a explicação de $X$ e $T$.

2. $R^2_{T \sim Z|X}$ é o R-quadrado parcial de $Z$ com relação à variável dependente $T$, o tratamento. Ou seja, quanto a variável omitida explica do tratamento, após descontar a explicação na variância dos controles.

3. $1 - R^2_{T \sim Z|X}$ é o complemento do R-quadrado parcial de Z sobre T, após controlar para $X$. Representa a variância não-explicada em $T$ por $Z$ após descontar a influência de $X$.

E os demais termos são os desvios-padrões dos resíduos da regressão restrita, sem a variável não-observada $Z$ e da regressão do tratamento em função dos controles. Em nosso exemplo, nós observamos a variável regilião evangélica, mas vamos supor que não é calcular o viés de variável omitida pela fórmula tradicional.

```{r vies r-squared, echo=TRUE, message=FALSE}

###
# tratamento é ideologia.
# omitida é religião

beta2 <- coef(reg_full)[6]
cov_x1x2 <- cov(br_latbar_23$ideologia, br_latbar_23$evangelico)
var_x1 <- var(br_latbar_23$ideologia)
vies <- beta2*cov_x1x2/var_x1
print(vies)

coef_res_sem_vies = coef(reg_res_evan)[2] - vies
all.equal(coef_res_sem_vies, coef(reg_full)[2])
coef_res_sem_vies - coef(reg_full)[2]

```
O que essa equação nos mostra é que o viés depende de quanto a variável omitida explica o tratamento e quanto o tratamento é explicado pela variável omitida. O que é bastante intuitivo. E para entender melhor essa equação, vale falar o viés de relativo, definido como a razão entre o viés estimado e o coeifciente do tratamento estimado.

$$
Viés relatrivo = \frac{\hat{<viés>}}{\hat{\beta_1}} = \frac{f(R^2_{y\sim Z|T,X},  R^2_{D \sim Z|X)})}{t_-do_-tratamento \cdot \sqrt(df)}
$$

a parcela da variável resposta que o tratamento explica (medida pelo R-quadrado parcial do tratamento) dá uma medida da robustez da regressão à potenciais vieses de variável omitida. Uma intuição simples é que se o tratamento explicar 100% da variação, então não há, por definição, viés de variável omitida. Se explicar 95%, há muito pouco que um viés de variável omitida pode viesar as estimativas.

```{r analise sensibilidade, echo=TRUE, message=FALSE}
library(sensemakr)

# loads dataset
data("darfur")

# runs regression model
model <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
                         pastvoted + hhsize_darfur + female + village, data = darfur)

summary(model)
# runs sensemakr for sensitivity analysis
sensitivity <- sensemakr(model = model, 
                         treatment = "directlyharmed",
                         benchmark_covariates = "female",
                         kd = 1:3)

# short description of results
sensitivity

summary(sensitivity)
```

## Referências

https://grantmcdermott.com/interaction-effects/

Brambor, T., Clark, W. R., & Golder, M. (2006). Understanding interaction models: Improving empirical analyses. Political analysis, 14(1), 63-82.

Corrêa, D. S. (2015). Economy, Ideology and Elections in Latin America. DADOS: Revista de Ciencias Sociais, 58(2), 401.