• Pessoas que trabalham com análise de dados frequentemente são demandadas para realizar previsões, não é tarefa fácil, mas alguém tem que fazer. Nessa unidade vamos falar sobre previsão.
• Mais uma vez usaremos de exemplo para introduzir o assunto, no caso previsões nas eleições americanas.
• Nessa unidade trataremos de regressões lineares, uma ferramenta essencial para cientistas sociais que trabalham com dados.

• O sistema de colégio eleitoral dos EUA, único no mundo, faz com que prever as eleições americanas seja particularmente desafiador. Um candidato ganha quando consegue a maioria no colégio aleitoral, não do voto popular.
• Em 2008 o colégio eleitoral tinha 538 eleitores, desses 535 pertenciam aos 50 estados e os três restantes ao Distrito de Columbia (Washington DC). Na maioria dos casos o eleitor do colégio eleitoral vota no candidato que ganhou a maioria dos votos no estado que ele representa,em alguns estados não fazer isso é crime.
• Desta forma as eleições ocorrem em sistema de o ganhador leva tudo, quem tem a maioria no estado leva todos os votos desse estado no colégio eleitoral. Em 2008, Obama teve 365 votos no colégio eleitoral e McCain ficou com 173 votos.

• Por conta do sistema de colégio eleitoral, prever as eleições nos EUA equivale a prever as eleições em cada estado. Uma vitória por pequena margem em um estado pode definir as eleições, de fato isso aconteceu nas eleições de 2000 quando Bush venceu Al Gore após a recontagem de votos na Flórida.
• Antes de começar com o exemplo, será útil tratar de loops e discutir melhor o uso de condições no R.

• Em muitos casos desejamos repetir uma operação várias vezes com pequenas alterações em cada repetição. Por exemplo, para fazer previsão nas eleições americanas será necessário repetir os procedimentos para cada um dos 50 estados. Não é razoável escrever o mesmo código 50 vezes, ainda que você tope a empreitada imagine um estudo que envolva os 5.570 municípios do Brasil.

• Um loop é uma estrutura de programação que permite repetir operações semelhantes com um código compacto.
• No R o loop do tipo for (i in X) vai criar um loop onde \(i\) é o contador (counter) e \(X\) é o vetor de valores que \(i\) irá tomar.

• A estrutura de um loop do tipo for no **R* é:

for(i in X){
  expressão 1
  expressão 2
  ...
  expressão N
}

• As expressões de 1 a N serão repetidas para cada valor \(i\) do vetor \(X\).
• Em cada iteração, \(i\) terá o valor correspondente no vetor \(X\), começando do primeiro elemento e terminando no último.
• Façamos um exemplo para multiplicar cada elemento de um vetor por 2.

• É boa prática criar um vetor “vazio” onde serão armazenados os resultados da operação. Vamos fazer isso usando a função rep() para criar um vetor de NAs.

values <- c(2,4,6)
values

## [1] 2 4 6

n <- length(values)
n

## [1] 3

results <- rep(NA,n)
results

## [1] NA NA NA

• O loop deve pegar cada elemento do vetor values, multiplicar por 2 e escrever uma mensagem com o resultado. Para escrever o resultado vamos usar a função str_c() do pacote stringr. Essa função combina múltiplos objetos em um único vetor de caracteres. A função str_c() é semelhante à função paste0() do R básico, mas é mais compatível com a estrutura do tidyverse.
• Também será usada a função seq_along() para enumerar as iterações.

library(stringr)
x1 <- "Roberto"
x2 <- "Ellery"
str_c(x1,x2, sep = " ")

## [1] "Roberto Ellery"

x1 <- "ECO"
x2 <- "UnB"
str_c(x1,x2, sep = "/")

## [1] "ECO/UnB"

seq_along(values)

## [1] 1 2 3

seq_along(c("a", "b", "c", "d"))

## [1] 1 2 3 4

• Com essas funções podemos criar nosso loop:

for (i in seq_along(values)) {
  results[i] = values[i] * 2
  print(str_c(values[i], " vezes 2 é igual a ", results[i]))
}

## [1] "2 vezes 2 é igual a 4"
## [1] "4 vezes 2 é igual a 8"
## [1] "6 vezes 2 é igual a 12"

results

## [1]  4  8 12

• O loop do exemplo é desnecessário, pois poderíamos obter o mesmo resultado multiplicando o vetor values por dois:

values * 2

## [1]  4  8 12

• Sempre que possível evite usar loops, apesar de intuitivos eles são computacionalmente intensivos.

• Uma boa prática quando usando loop é usar a função print() para avisar o que está acontecendo em cada iteração, isso ajuda a identificar onde ocorre um erro. Para ilustrar, primeiro criemos um data.frame com os dados fictícios.

data <- data.frame("a" = 1:2, "b" = c("hi", "hey"), "c"=3:4)
data

##   a   b c
## 1 1  hi 3
## 2 2 hey 4

• Agora vamos fazer um loop com um erro na segunda iteração.

results <- rep(NA, 3)

for (i in seq_along(data)) {
  print(str_c("iteration ", i))
  results[i] <- median(data[,i])
}

results

results <- rep(NA, 3)

for (i in seq_along(data)) {
  print(str_c("iteration ", i))
  results[i] <- median(data[,i])
}

## [1] "iteration 1"
## [1] "iteration 2"

## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
## argumento não é numérico nem lógico: retornando NA

## [1] "iteration 3"

results

## [1] 1.5  NA 3.5

• Nas unidades anteriores falamos das funções if_else() e case_when(), ambas executam alguma tarefa quando uma condição é verdadeira. Nesta unidade vamos tratar dos comandos if () {} e if () {} else {}, esses comandos permitem condicionar um bloco de código a uma expressão ser verdadeira.

• A construção if () {} executa o código em {} quando a expressão em () é verdadeira, tem a forma:

if (X) {
  expressão 1
  expressão 2
  ...
  expressão N
}

• Se X for verdadeiro, as expressões de 1 a N serão executadas. Se X for falso, as expressões de 1 a N serão ignoradas.

• Exemplo

operation <- "add"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
}

if (operation == "multiply") {
  print("fazer multiplicação 4 * 4")
  4 * 4
}

operation <- "add"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
}

## [1] "fazer adição 4 + 4"

## [1] 8

if (operation == "multiply") {
  print("fazer multiplicação 4 * 4")
  4 * 4
}

• A construção if () {} else {} executa o código no primeiro {} quando a expressão em () é verdadeira, se a expressão em () for falsa o código no segundo {} será executado:

if (X) {
  expressão 1a
  expressão 2a
  ...
  expressão Na
} else {
  expressão 1b
  expressão 2b
  ...
  expressão Nb
}

• Se X for verdadeiro, as expressões de 1a a Na serão executadas. Se X for falso, as expressões de 1b a Nb serão executadas.

• Exemplo

operation <- "multiply"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
} else {
  print("fazer multiplicação 4 * 4")
  4 * 4
}

operation <- "multiply"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
} else {
  print("fazer multiplicação 4 * 4")
  4 * 4
}

## [1] "fazer multiplicação 4 * 4"

## [1] 16

• Podemos fazer construções condicionais ainda mais complexas usando else if ()

if (X) {
  expressão 1a
  ...
  expressão Na
} else if (Y) {
  expressão 1b
  ...
  expressão Nb
} else {
  expressão 1c
  ...
  expressão Nc
}

• Se X for verdadeiro, as expressões 1a a Na serão executadas, de Y for verdadeiro, as expressões 1b a Nb serão executadas, finalmente, se X e Y forem falsos, as expressões 1c a Nc serão executadas. Repare que a ordem das expressões importa.

• Exemplo

operation <- "subtract"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
} else if (operation == "multiply") {
  print("fazer multiplicação 4 * 4")
  4 * 4
} else {
  print(str_c("`", operation, "´ é inválido. Use `add´ ou `multiply´."))
}

operation <- "subtract"

if (operation == "add") {
  print("fazer adição 4 + 4")
  4 + 4
} else if (operation == "multiply") {
  print("fazer multiplicação 4 * 4")
  4 * 4
} else {
  print(str_c("`", operation, "´ é inválido. Use `add´ ou `multiply´."))
}

## [1] "`subtract´ é inválido. Use `add´ ou `multiply´."

• Os comandos condicionais podem ser usados dentro de um loop do tipo for. O exemplo abaixo avalia números inteiros de 1 a 5, se o número for par retorna o número mais o número (dobro do número), se for ímpar retorna o número multiplicado pelo número (quadrado).

values <- 1:5
n <- length(values)
results <- rep(NA, n)

for (i in seq_along(values)) {
  x <- values[i]
  r <- x %% 2
  if (r==0) {
    print(str_c(x, " é par, operação é adição: ", x, "+", x))
    results[i] <- x + x
  } else {
    print(str_c(x, " é ímpar, operação é multiplicação: ", x, "*", x))
    results[i] <- x * x
  }
}
results

## [1] "1 é ímpar, operação é multiplicação: 1*1"
## [1] "2 é par, operação é adição: 2+2"
## [1] "3 é ímpar, operação é multiplicação: 3*3"
## [1] "4 é par, operação é adição: 4+4"
## [1] "5 é ímpar, operação é multiplicação: 5*5"

## [1]  1  4  9  8 25

• Conhecendo melhor o loop e a execução de expressões condicionais podemos seguir com o exemplo de previsão das eleições nos EUA.
• Os dados em pres08 são dos resultados em cada estado e em polls08 estão resultados de pesquisas nos estados. Queremos usar os dados das pesquisas para prever o resultado das eleições.
• Depois de carregar os dados a primeira taerfa será criar em cada data.frame uma variável chamada margin com as diferenças entre Obama e McCain.

• Variáveis em pres08:
  • state: sigla do estado
  • state.name: nome do estado
  • Obama: percentual de votos para Obama
  • McCain: percentual de votos para McCain
  • EV: número de votos do estado no colégio eleitoral

• Variáveis em polls08:
  • state: sigla do estado
  • Obama: previsão do percentual de votos para Obama
  • McCain: previsão percentual de votos para McCain
  • pollster: nome da organização que fez a pesquisa
  • middate: período que a pesquisa foi feita (meia data)

library(tidyverse)

pres08 <- read.csv("pres08.csv")
polls08 <- read.csv("polls08.csv")

pres08 <- pres08 %>%
  mutate(margin = Obama - McCain)

polls08 <- polls08 %>%
  mutate(margin = Obama - McCain)

head(pres08)

##   state.name state Obama McCain EV margin
## 1    Alabama    AL    39     60  9    -21
## 2     Alaska    AK    38     59  3    -21
## 3    Arizona    AZ    45     54 10     -9
## 4   Arkansas    AR    39     59  6    -20
## 5 California    CA    61     37 55     24
## 6   Colorado    CO    54     45  9      9

head(polls08)

##   state         Pollster Obama McCain    middate margin
## 1    AL      SurveyUSA-2    36     61 2008-10-27    -25
## 2    AL Capital Survey-2    34     54 2008-10-15    -20
## 3    AL      SurveyUSA-2    35     62 2008-10-08    -27
## 4    AL Capital Survey-2    35     55 2008-10-06    -20
## 5    AL      Rasmussen-1    39     60 2008-09-22    -21
## 6    AL      SurveyUSA-2    34     64 2008-09-16    -30

• Para cada estado será criada uma variável com margem para o Obama prevista nas pesquisas mais próximas da eleição, ou seja, vamos calcular a média das pesquisas feitas em cada estado no dia mais próximo da eleição.
• O dia da última pesquisa pode ser diferente entre os estados e também é possível que exista mais de uma pesquisa feita no mesmo dia (para ser preciso com a mesma meia data) em um determinado estado.
• Para fazer isso vamos criar um vetor vazio de tamanho 51 (os 50 estados e DC) que será chamado poll.pred, depois faremos um loop onde cada iteração tratará dos dados de um estado.

• Antes, porém, será preciso falar um pouco sobre como o R trabalha com datas. No R existem classes específicas para datas, a mais básica é chamada Date.

• O pacote lubridate contém funções úteis para trabalhar com datas, entre elas a função ymd() que transforma caracteres em datas com formato year-month-day. Existem outras funções semelhantes para diferentes padrões de data, por exemplo ymd() e dmy().

library(lubridate)

x <- ymd("2008-11-04")
y <- ymd("2008/9/1")

class(x)

## [1] "Date"

class(y)

## [1] "Date"

x

## [1] "2008-11-04"

y

## [1] "2008-09-01"

subtraction <- x - y
subtraction

## Time difference of 64 days

class(subtraction)

## [1] "difftime"

as.numeric(subtraction)

## [1] 64

• Com essas ferramentas vamos criar uma variável com os dias entre a pesquisa e as eleições, o nome da variável será DaysToElection.
• Na sequência vamos calcular as médias dos resultados das pesquisas mais próximas das eleições.
• A função unique() será usada para obter a sigla única de cada estado.

• No data.frame polls08 a variável middate não é da classe Date

str(polls08)

## 'data.frame':    1332 obs. of  6 variables:
##  $ state   : chr  "AL" "AL" "AL" "AL" ...
##  $ Pollster: chr  "SurveyUSA-2" "Capital Survey-2" "SurveyUSA-2" "Capital Survey-2" ...
##  $ Obama   : int  36 34 35 35 39 34 36 25 35 34 ...
##  $ McCain  : int  61 54 62 55 60 64 58 52 55 47 ...
##  $ middate : chr  "2008-10-27" "2008-10-15" "2008-10-08" "2008-10-06" ...
##  $ margin  : int  -25 -20 -27 -20 -21 -30 -22 -27 -20 -13 ...

• Vamos usar a função ymd() para resolver esse problema.

polls08 <- polls08 %>%
  mutate(middate = ymd(middate))

str(polls08)

## 'data.frame':    1332 obs. of  6 variables:
##  $ state   : chr  "AL" "AL" "AL" "AL" ...
##  $ Pollster: chr  "SurveyUSA-2" "Capital Survey-2" "SurveyUSA-2" "Capital Survey-2" ...
##  $ Obama   : int  36 34 35 35 39 34 36 25 35 34 ...
##  $ McCain  : int  61 54 62 55 60 64 58 52 55 47 ...
##  $ middate : Date, format: "2008-10-27" "2008-10-15" ...
##  $ margin  : int  -25 -20 -27 -20 -21 -30 -22 -27 -20 -13 ...

election_date <- ymd("2008-11-04")

polls08 <- polls08 %>%
  mutate(DaysToElection = as.numeric(election_date - middate))

poll.pred <- rep(NA,51)

st.names <- unique(polls08$state)

names(poll.pred) <- as.character(st.names)

for(i in seq_along(st.names)) {
  state.data <- polls08 %>%
    filter(state == st.names[i])
  
  min_days <- min(state.data$DaysToElection)
  
  state.data <- state.data %>%
    filter(DaysToElection == min_days)
  
  poll.pred[i] <- mean(state.data$margin)
}

head(poll.pred, n=10)

##    AL    AK    AZ    AR    CA    CO    CT    DC    DE    FL 
## -25.0 -19.0  -2.5  -7.0  24.0   7.0  25.0  69.0  30.0   2.0

• Quão acurada é nossa previsão? Para responder essa pergunta vamos comparar o previsto com o observado em cada estado.
• A primeira medida de erro que vamos usar é chamada de erro de predição e consiste na diferença entre o que aconteceu e o que foi previsto.
• Após calcular o erro de previsão em cada estado vamos calcular a média desses erros.

• É importante ter certeza que os estados estão na mesma ordem, mais a frente vamos aprender a juntar bases de dados o que resolve problemas desse tipo.

errors <- pres08$margin - poll.pred
head(errors, n=10)

##    AL    AK    AZ    AR    CA    CO    CT    DC    DE    FL 
##   4.0  -2.0  -6.5 -13.0   0.0   2.0  -2.0  16.0  -5.0   1.0

mean(errors)

## [1] 1.062092

• Um problema de trabalhar com médias das diferenças entre o que ocorreu e o que foi previsto é que erros para mais e erros para menos podem ser cancelar, isso pode passar uma falsa impressão de acuracia.
• Para não correr esse risco podemos usar a raiz da média dos quadrados do erro, como o quadrado de um número negativo é positivo não ocorre de erros para baixo cancelarem erros para cima ou vice-versa.

sqrt(mean(errors^2))

## [1] 5.90894

• O resultado indica que, em média, a magnitude do erro das pesquisas foi de aproximadamente 6%.

• O erro de previsão é definido como: \[\mbox{erro de predição}=\mbox{valor observado}-\mbox{valor previsto}\]
• A média dos erros de predição é chamada viés, uma previsão é dita não-viesada quando o viés é zero. A raiz da média do quadrado dos erros é chamada de raiz do erro quadrado médio e representa a média da magnitude do erro de previsão.

• Podemos analisar toda a distribuição do erro usando um histograma.

errors_tib <- as_tibble(errors)

errors_tib %>%
  ggplot(aes(errors, after_stat(density))) +
  geom_histogram(binwidth = 5, boundary = 0) +
  geom_vline(xintercept = mean(errors_tib$value)) +
  annotate("text", x=-.2, y=0.07, hjust=1, label = "average error", size=14/.pt) +
  labs(title = "Poll prediction error", y = "Density",
       x= "Error in predicit margin for Obama (percentage points") +
  theme_classic(base_size = 13)

• O histograma mostra uma grande variação do erro entre os estados, porém na maioria dos casos o erro é relativamente pequeno, erros grandes são menos frequentes. A distribuição tem um formato de sino em volta de zero.

• Para examinar com mais cuidado o erro em cada estado podemos fazer um gráfico de dispersão com a previsão para o estado no eixo horizontal e o resultado do estado no eixo vertical.
• Um estado onde a previsão é igual ao resultado ficará na linha de 45o (que estará no gráfico). Quando a previsão errar a favor do Obama (margem prevista maior do que observada) o estado fica acima da linha de 45o, quando a previsão errar a favor do McCain o estado fica abaixo da linha.
• Para fazer o gráfico será útil colocar os resultados previstos e observados na mesma base de dados, faremo sisso com a função cbind(). Para usar cbind() as observações devem estar na mesma ordem nas colunas que serão combinadas.

pres08 <- pres08 %>%
  cbind(poll.pred = poll.pred)

ggplot(data=pres08, aes(x=poll.pred, y=margin)) +
  geom_text(aes(label = state)) +
  geom_abline(linetype="dashed") +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  ylim(-40, 100) +
  xlim(-40, 100) +
  labs(x = "Poll results", y="Acutal election results")

• Washington DC e Vermont, DC e VT no gráfico, estão bem acima da linha de 45o, ou seja, as pesquisas erraram a favor de Obama. No Arizona e na Luisiana, AR e LA no gráfico, aconteceu o contrário, ambos estão bem abaixo da linha mostrando que as pesquisas erraram a favor do McCain.
• Dado o modelo de eleiçõs nos EUA, erros desse tipo não impactam a previsão sobre o número de votos no colégio eleitoral. Por outro lado, em um estado muito disputado, é possível que um erro pequeno atrapalhe a previsão de votos no colégio eleitoral. Suponha que a previsão foi de vitória de McCain por menos de 1% e o ocorrido foi vitória de Obama por menos de 1%, na grande maioria dos estados a previsão terá errado todos os votos do estado no colégio eleitoral.

• Existem dois tipos de erros de predição quando as previsões apontam o vencedor errado.
• No nosso exemplo, no quadrante superior esquerdo a previsão é que Obama perderia (a esquerda de zero/linha vertical), mas Obama ganhou (acima de zero na linha/horizontal). No quadrante inferior direito ocorre o contrário, a previsão é de vitória de Obama (a direita de zero/linha vertical), mas Obama perdeu (abaixo de zero/linha horizontal).
• No gráfico fica claro que na grande maioria dos estados as pesquisas acertaram o vencedor (ainda que tenham errado na margem), mas em três as pesquisas erraram. Nesses três estados a disputa foi tão apertada que é difícil ver no gráfico o que ocorreu.

• A função sign() retorna o sinal de um número, especificamente retorna -1 se o número for negativo e 1 se o número for positivo. Quando a margem tem mesmo sinal na pesquisa e no resultado observado é porque a pesquisa acertou o vencedor (no gráfico correponde aos quadrante inferior esquerdo e superior direito). Quando os sinais são diferentes a pesquisa errou o vencedor.
• Vamos usar a função sign() para encontrar os estados onde as pesquisas erraram.

pres08 <- pres08 %>%
  mutate(correct = if_else(sign(poll.pred) == sign(margin), 1, 0))

pres08 %>% filter(correct == 0) %>%
  select(state.name, Obama, McCain, margin, poll.pred, correct)

##        state.name Obama McCain margin poll.pred correct
## IN        Indiana    50     49      1        -5       0
## MO       Missouri    48     49     -1         1       0
## NC North Carolina    50     49      1        -1       0

• O problema de prever o resultado por categoria ou classe é chamado de classificação. No nosso exemplo, para cada estado queremos prever se Obama ganha ou não.
• Em problemas de classificação a predição é claramente certa ou errada. Uma previsão errada é chamada de erro de classificação (missclassification).
• No exemplo, a taxa de erro de classificação foi de 3/51 (3 erros em 51 previsões), aproximadamente 6%.

• Em um problema de classficação binário existem dois tipos de erros de classificação. Podemos prever que Obama vence em um estado e ele perde ou podemos prever que Obama perde e ele vence.
• Se chamarmos a vitória de Obama de um positivo, então o primeiro tipo de erro é um falso positivo (erramos ao prever que ele ganharia, ou seja, a previsão era positivo e deu negativo) e o segundo tipo de erro é um falso negativo (a previsão era negativo e deu positivo).

• No exemplo, Missouri foi um falso positivo enquanto Indiana e Carolina do Norte foram falsos negativos.
• A matriz de confusão (confusion matrix) resume os erros e acertos na classificação.

• Classificação é o problema de prever um resultado por categorias. Na classficação a previsão é correta ou incorreta. Em um problema de classficação binária existem dois tipos de erros de classficação: falsos positivos e falsos negativo, esses erros representam erroneamente prever o resultado como positivo ou negativo, respectivamente.

• Podemos calcular o número de votos para Obama no colégio eleitoral previsto pelas pesquisas e comparar com o número de votos que ele efetivamente obteve. Para ser eleito Obama precisava de 270 votos.

#Votos obtidos
pres08 %>%
  filter(margin > 0) %>%
  summarize(total_EV = sum(EV))

##   total_EV
## 1      364

#Votos previstos
pres08 %>%
  filter(poll.pred > 0) %>%
  summarize(pred_EV = sum(EV))

##   pred_EV
## 1     349

• Nos EUA o voto popular não define o vencedor das eleições, mas pode ser visto como uma medida da opinião pública em relação aos candidatos.
• Com essa perspectiva, podemos usar as pesquisas para avaliar como a opinião pública se comportou nos meses anteriores às eleições.
• Os dados estão em pollsUS08, as variáveis são as mesmas que estamos usando para as pesquisas. Para cada um dos 90 dias anteriores às eleições calcularemos as médias de cada candidato em todas as pesquisas realizadas na semana anterior ao dia. A eleição foi no dia 4/11/2008.

pollsUS08 <- read.csv("pollsUS08.csv")
str(pollsUS08)

## 'data.frame':    526 obs. of  4 variables:
##  $ Pollster: chr  "IBD/TIPP " "GWU (Lake/Tarrance) " "Diageo/Hotline " "Newsweek " ...
##  $ McCain  : int  48 51 43 44 44 45 49 42 44 48 ...
##  $ Obama   : int  36 39 38 46 47 47 42 48 44 48 ...
##  $ middate : chr  "2007-07-04" "2007-07-11" "2007-07-14" "2007-07-19" ...

head(pollsUS08, n=15)

##                Pollster McCain Obama    middate
## 1             IBD/TIPP      48    36 2007-07-04
## 2  GWU (Lake/Tarrance)      51    39 2007-07-11
## 3       Diageo/Hotline      43    38 2007-07-14
## 4             Newsweek      44    46 2007-07-19
## 5            Rasmussen      44    47 2007-07-19
## 6             ABC/Post      45    47 2007-07-19
## 7                 Time      49    42 2007-07-24
## 8             Newsweek      42    48 2007-07-26
## 9            Rasmussen      44    44 2007-08-09
## 10    USA Today/Gallup      48    48 2007-08-11
## 11         WNBC/Marist      44    41 2007-08-15
## 12          Quinnipiac      43    43 2007-08-17
## 13               Zogby      40    44 2007-08-24
## 14                Time      42    46 2007-08-26
## 15            Newsweek      43    45 2007-08-30

pollsUS08 <- pollsUS08 %>%
  mutate(middate = ymd(middate))

str(pollsUS08)

## 'data.frame':    526 obs. of  4 variables:
##  $ Pollster: chr  "IBD/TIPP " "GWU (Lake/Tarrance) " "Diageo/Hotline " "Newsweek " ...
##  $ McCain  : int  48 51 43 44 44 45 49 42 44 48 ...
##  $ Obama   : int  36 39 38 46 47 47 42 48 44 48 ...
##  $ middate : Date, format: "2007-07-04" "2007-07-11" ...

all_dates <- seq(min(pollsUS08$middate), election_date, by="days")
head(all_dates)

## [1] "2007-07-04" "2007-07-05" "2007-07-06" "2007-07-07" "2007-07-08"
## [6] "2007-07-09"

tail(all_dates)

## [1] "2008-10-30" "2008-10-31" "2008-11-01" "2008-11-02" "2008-11-03"
## [6] "2008-11-04"

prior_days <- 7

vote_avg <- vector(length(all_dates), mode = "list")

for (i in seq_along(all_dates)) {
  date <- all_dates[i]
  week_data <-
    filter(pollsUS08, as.numeric(date - middate) >= 0, 
           as.numeric(date - middate) < prior_days) %>%
    summarize(Obama = mean(Obama, na.rm = TRUE),
              McCain = mean(McCain, na.rm = TRUE))
  week_data$date <- date
  vote_avg[[i]] <- week_data
}

head(vote_avg, n=4)

## [[1]]
##   Obama McCain       date
## 1    36     48 2007-07-04
## 
## [[2]]
##   Obama McCain       date
## 1    36     48 2007-07-05
## 
## [[3]]
##   Obama McCain       date
## 1    36     48 2007-07-06
## 
## [[4]]
##   Obama McCain       date
## 1    36     48 2007-07-07

• Para facilitar nossa vida vamos tranformar vote_avg de um objeto da classe list para um objeto da classe data.frame.

vote_avg_df <- bind_rows(vote_avg)
head(vote_avg_df, n=4)

##   Obama McCain       date
## 1    36     48 2007-07-04
## 2    36     48 2007-07-05
## 3    36     48 2007-07-06
## 4    36     48 2007-07-07

• O próximo passo é fazer uma figura com os resultados dos 90 dias anteriores à eleição.

vote_avg_df %>%
  filter(election_date - date <= 90) %>%
  ggplot() +
  geom_point(aes(x=date, y=Obama), color = "blue", shape = 1) +
  geom_point(aes(x=date, y=McCain), color = "red", shape=1) +
  ylim(40,55) +
  labs(y="Support for candidate (percentage points)", x="Date") +
  annotate("text", x=ymd("2008-08-15"), y=47, label="Obama", color="blue") +
  annotate("text", x=ymd("2008-08-15"), y=41, label="McCain", color="red") +
  geom_vline(xintercept = election_date) +
  geom_point(aes(x=election_date, y=52.93), color="blue") +
  geom_point(aes(x=election_date, y=45.65), color="red") +
  theme_classic()

• Até agora fizemos previsões tomando médias dos resultados de pesquisas eleitorais. Uma alternativa é fazer previsões usando modelos estatísticos. O modelo de regressão linear é um dos mais básicos modelos estatísticos e não pode faltar na caixa de ferramentas de cientistas sociais aplicados.

• Estudos de psicologia apontam que prever o resultado de eleições com base no rosto do candidato é melhor do que apostar no acaso.
• A pesquisa consiste em mostrar retratos em preto e branco de dois candidatos que disputam um cargo (nos EUA as elições para deputado seguem o modelo distrital, é comum que sejam disputados entre dois candidatos) por menos de um segundo e pedir para o entrevistado avaliar os dois candidatos com base na competência.
• Os entrevistados não conhecem os candidatos, em particular não sabem o partido do candidato nem se o candidato está disputando reeleição.

• Os pesquisadores usam a medida de competência para prever o resultado da eleição. A medida de competência é proporção de entrevistados que classificam o candidato como mais competente do que o outro.
• A hipótese testada é se uma avaliação rápida da aparência pode prever resultados eleitorais.
• O exemplo tem por base: Alexander Todorov, Anesu Mandisodza, Amir Goren e Crystal Hall. Inference of competence from faces predict election outcome; Science, v.308, n.10, jun, 2005.
• Os dados estão do data.frame face.

• Variáveis:
  • congress: número da legislatura
  • year: ano da eleição
  • state: estado da eleição
  • winner: vencedor
  • loser: perdedor
  • w.party: partido do vencedor
  • l.party: partido do perdedor
  • d.votes: número de votos para o candidato Democrata
  • r.votes: número de votos para o candidato Republicano
  • d.comp: medida de competência do candidato Democrata
  • r.comp: medida de competência do candidato Republicano

face <- read.csv("face.csv")

• Para realizar a análise os autores criaram uma variável para capturar a margem do candidato democrata. Para o candidato de cada partido foi calculada a proporção de votos do candidato no total de votos para os candidatos dos dois partidos.
• A margem é a proporção do candidato Democrata menos a proporção do candidato Republicano. Se a margem for positiva indica vitória do Democrata, se for negativa a vitória foi do Republicano.
• O primeiro exercício consiste em um gráfico de dispersão com a margem do eixo vertical e a medida de competência do candidato democrata no eixo horizontal.

• Criar a variável com a margem.

face <- mutate(face,
               d.share = d.votes/(d.votes + r.votes),
               r.share = r.votes/(d.votes + r.votes),
               diff.share = d.share - r.share)

face %>%
  ggplot(aes(d.comp, diff.share, color=w.party)) +
  geom_point() +
  scale_colour_manual(values = c(D = "blue", R = "red")) +
  labs(x="Competence score for Democrats",
       y="Democrats margin in vote share",
       title = "Facial competence and vote share") +
  ylim(-1,1) +
  xlim(0,1) +
  theme(legend.position = "none")

• Na unidade anterior vimos que a correlação representa o grau de associação entre duas variáveis. Uma correlação positiva signfica que uma variável é mais provável de estar acima da média quando a outra variável está acima da média, vale o contrário para correlação negativa.
• A “nuvem” de dados com inclinação positiva na figura anterior indica uma correlação positivia entre a medida de competência e margem de vitória do candidato Democrata.

cor(face$d.comp, face$diff.share)

## [1] 0.4327743

• Correlação mede relação linear entre duas variáveis, uma “nuvem” de dados com inclinação positiva indica uma correlação positiva, uma nuvem com inclinação negativa indica uma correlação negativa.
• Quanto mais próxima de 1 a correlação, mais visível será a inclinação positiva da nuvem, vale o mesmo para correlações negativas.
• A próxima figura ilustra esse ponto.

• Repare que no gráfico abaixo e a esquerda a correlação é próxima de zero, porém existe uma relação entre as variáveis. Ocorre que não é uma relação linear. O exemplo mostra que ausência de correlação não implica em ausência de relação entre as variáveis.
• Análise de dados tem suas armadilhas, o pacote datasauRus contém bases de dados onde a relação entre duas variáveis \(x\) e \(y\) apresentam média, desvio padrão e correlação entre \(x\) e \(y\) muito semelhantes, mas…

library(datasauRus)

data("datasaurus_dozen")

datasaurus_dozen %>%
  filter(dataset %in% c("bullseye", "dino", "circle", "x_shape")) %>%
  group_by(dataset) %>%
  summarise(media_x = mean(x),
            media_y = mean(y),
            sd_x = sd(x),
            sd_y = sd(y),
            cor_xy = cor(x,y))

## # A tibble: 4 × 6
##   dataset  media_x media_y  sd_x  sd_y  cor_xy
##   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>
## 1 bullseye    54.3    47.8  16.8  26.9 -0.0686
## 2 circle      54.3    47.8  16.8  26.9 -0.0683
## 3 dino        54.3    47.8  16.8  26.9 -0.0645
## 4 x_shape     54.3    47.8  16.8  26.9 -0.0656

• O coeficiente de correlação quantifica a relação linear entre duas variáveis. A “nuvem” de dados no gráfico de dispersão apresentando uma inclinação positiva é sinal de correlação positiva, se a inclinação for negativa é sinal de correlação negativa. Via de regra, correlação não é adequada para representar relações não lineares.

• A correlação é uma medida de relação linear entre duas variáveis, mas muitas vezes queremos carcaterizar melhor essa relação. Para isso usamos um modelo linear: \[Y=\alpha+\beta X+\varepsilon\]
• \(Y\) é a variável de resultado, também chamada de resposta ou dependente, \(X\) é a variável preditora, também chamada de explicativa ou independente. No nosso exemplo, a medida de competência é a variável preditora, \(X\), e a diferença nas frações de votos dos candidatos é a variável de resultado, \(Y\).

• Toda reta pode ser definida pelo intercepto, \(\alpha\), e a inclinação, \(\beta\). O intercepto representa o valor médio de \(Y\) quando \(X\) é zero, a inclinação representa o aumento médio de \(Y\) quando \(X\) aumenta de uma unidade. O intercepto e a inclinação são chamados de coeficientes.
• O termo de erro, \(\varepsilon\), captura desvios de uma relação perfeitamente linear entre as variáveis \(Y\) e \(X\).

• Modelos como o linear são usados para aproximar o processo gerador de dados.
• O livro cita uma máxima do estatístico George Box de que “todos os modelos são errados, mas alguns são úteis”. Esse é um ponto fundamental também em economia, melhor do que buscar por modelos verdadeiros, que pode se tornar uma busca pelo Santo Graal, é buscar por modelo úteis. Lembrem do as if do Friedman.

• Como não conhecemos os valores de \(\alpha\) e \(\beta\), usamos os dados para estimar esses valores. Por convenção, os valores estimados são denotados pela letra do parâmetro desconhecido com um chapeu, no nosso modelo, o valor estimado para \(\alpha\) será denotado por \(\hat{\alpha}\) e o valor estimado para \(\beta\) será denotado por \(\hat{\beta}\).

• Uma vez conhecidos \(\hat{\alpha}\) e \(\hat{\beta}\), podemos obter a reta de regressão usando \(\hat{\alpha}\) como intercepto e \(\hat{\beta}\) como inclinação. Para um dado valor de \(X\), chame de \(x\) a reta de regressão permite obter o valor previsto, ou valor ajustado (fitted value) de \(Y\): \[\hat{Y}=\hat{\alpha}+\hat{\beta}x\]
• A diferença entre o valor observado da variável dependente, \(Y\), e o valor previsto desta variável \(\hat{Y}\) é chamado de resíduo ou erro de previsão e é dado por: \[\hat{\varepsilon}=Y-\hat{Y}\]

• O modelo de regressão linear é defindo por: \[Y=\alpha+\beta X + \varepsilon\] Onde \(Y\) é a variável de resultado (ou resposta), \(X\) é a variável preditora (ou independente, ou ainda explicativa), \(\varepsilon\) é o erro ou termo de disturbância e \((\alpha, \beta)\) são os coeficientes. O parâmetro de inclinação, \(\beta\) representa o aumento médio em \(Y\) associado ao aumento de uma unidade em \(x\). Uma vez que as estimativas dos coeficientes \((\hat{\alpha}, \hat{\beta})\) são obtidas dos dados, podemos prever \(Y\) usando um dado valor do preditor, chame esse valor de \(x\), o valor previsto de \(Y\) é dado por \(\hat{Y} = \hat{\alpha} + \hat{\beta}x\). A diferença entre o valor observado e o valor previsto (ou ajustado) de \(Y\) é chamada de resíduo e denotada por \(\hat{\varepsilon} = Y - \hat{Y}\).

• A função lm() do R faz regressão linear. Essa função tem como principal argumento uma fórmula, do tipo Y ~ X, também é comum usar o argumento data para indicar onde estão os dados. No R o símbolo ~ denota uma fórmula.
• Se os dados estiverem na área de trabalho não é preciso usar o argumento data para indicar a origem dos dados, da mesma forma é possível dizer onde estão os dados usando o nome do data.frame seguindo de $ e do nome da variável.

fit <- lm(diff.share ~ d.comp, data = face)
fit

## 
## Call:
## lm(formula = diff.share ~ d.comp, data = face)
## 
## Coefficients:
## (Intercept)       d.comp  
##     -0.3122       0.6604

lm(face$diff.share ~ face$d.comp)

## 
## Call:
## lm(formula = face$diff.share ~ face$d.comp)
## 
## Coefficients:
## (Intercept)  face$d.comp  
##     -0.3122       0.6604

• O valor estimado para o intercepto diz que, em média, se nenhum entrevistado apontar o candidato democrata como mais competente do que o Republicano, o candidato democrata perde com uma diferença de margem de 0,31.
• O valor estimado para inclinação, o coeficiente de d.comp, diz que se a competência percebida aumentar em 10 pontos a diferença da margem aumenta em aproximadamente 6,6%.

• O objeto fit é da classe lm e possui várias informações sobre o resultado da função lm()

class(fit)

## [1] "lm"

names(fit)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

• Podemos acessar cada uma dessas informações pelo nome.

coef(fit)

## (Intercept)      d.comp 
##  -0.3122259   0.6603815

fit$coef

## (Intercept)      d.comp 
##  -0.3122259   0.6603815

head(fitted(fit))

##           1           2           3           4           5           6 
##  0.06060411 -0.08643340  0.09217061  0.04539236  0.13698690 -0.10057206

• Acessar os resultados pelo nome é útil, mas pode ser pouco prático em trabalhos mais intensos.
• O pacote broom que está na coleção de pacotes tidymodels permite transformar as informações guardadas em objetos de classe **lm* (e de outras classes também) em data.frames.

• A função glance() retorna um data.frame com estatísiticas importantes do modelo (mais na frente falamos sobre essas estatísticas).

library(tidymodels)
glance(fit)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic     p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>       <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.187         0.180 0.266      27.0 0.000000885     1  -10.5  27.0  35.3
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

• A função tidy() retorna um data.frame onde cada linha traz a estimativa e estatísticas sobre um coeficiente.

tidy(fit)

## # A tibble: 2 × 5
##   term        estimate std.error statistic     p.value
##   <chr>          <dbl>     <dbl>     <dbl>       <dbl>
## 1 (Intercept)   -0.312    0.0660     -4.73 0.00000624 
## 2 d.comp         0.660    0.127       5.19 0.000000885

• A função augment() retorna os dados originais, os valores previstos (ajustados), os resíduos e outras estatísticas relacionadas as observações.

augment(fit)

## # A tibble: 119 × 8
##    diff.share d.comp .fitted  .resid    .hat .sigma  .cooksd .std.resid
##         <dbl>  <dbl>   <dbl>   <dbl>   <dbl>  <dbl>    <dbl>      <dbl>
##  1     0.210   0.565  0.0606  0.150  0.00996  0.267 0.00160       0.564
##  2     0.119   0.342 -0.0864  0.206  0.0129   0.267 0.00394       0.778
##  3     0.0499  0.612  0.0922 -0.0423 0.0123   0.268 0.000158     -0.160
##  4     0.197   0.542  0.0454  0.151  0.00922  0.267 0.00151       0.570
##  5     0.496   0.680  0.137   0.359  0.0174   0.266 0.0163        1.36 
##  6    -0.350   0.321 -0.101  -0.249  0.0143   0.267 0.00644      -0.941
##  7     0.697   0.404 -0.0456  0.743  0.00979  0.258 0.0388        2.80 
##  8     0.266   0.603  0.0860  0.180  0.0118   0.267 0.00276       0.681
##  9    -0.372   0.539  0.0434 -0.415  0.00914  0.265 0.0113       -1.57 
## 10     0.0106  0.869  0.262  -0.251  0.0426   0.267 0.0206       -0.963
## # ℹ 109 more rows

• Podemos colocar a reta de regressão no gráfico usando a função geom_abline().

ggplot() +
  geom_point(data=face, aes(d.comp, diff.share), shape=1) +
  geom_abline(slope = coef(fit)["d.comp"],
              intercept = coef(fit)["(Intercept)"]) +
  scale_x_continuous("Competence score for Democrats",
                     breaks = seq(0,1,by=0.2), limits = c(0,1)) +
  scale_y_continuous("Democratic margin in vote shares",
                     breaks = seq(-1,1,by=0.5), limits = c(-1,1)) +
  geom_vline(xintercept=mean(face$d.comp), linetype="dashed") +
  geom_hline(yintercept=mean(face$diff.share), linetype="dashed") +
  ggtitle("Facial competence and vote share") +
  theme_classic()

• Outra maneira de colocar a reta de regressão no gráfico é usar a função geom_smooth().

ggplot(data=face, aes(d.comp, diff.share)) +
  geom_point(shape=1) +
  geom_smooth(method="lm", se=FALSE, color="black") +
  scale_x_continuous("Competence score for Democrats",
                     breaks = seq(0,1,by=0.2), limits = c(0,1)) +
  scale_y_continuous("Democratic margin in vote shares",
                     breaks = seq(-1,1,by=0.5), limits = c(-1,1)) +
  geom_vline(xintercept=mean(face$d.comp), linetype="dashed") +
  geom_hline(yintercept=mean(face$diff.share), linetype="dashed") +
  ggtitle("Facial competence and vote share") +
  theme_classic()

• A reta de regressão é a reta com “melhor ajuste” porque minimiza a magnitude do erro de previsão. Um método comum para estimar os parâmetros (intercepto e inclinação) é o dos mínimos quadrados. A ideia é escolher \(\hat{\alpha}\) e \(\hat{\beta}\) de forma a minimizar a soma dos quadrados dos resíduos (SSR, do inglês sum of squared residual): \[\mbox{SSR}=\sum_{i=1}^n\hat{\varepsilon}^2=\sum_{i=1}^n(Y_i-\hat{Y}_i)^2=\sum_{i=1}^n(Y_i-\hat{\alpha}-\hat{\beta}X_i)^2\]

• O valor do SSR pode ser difícil de interpretar, mas é próximo do conceito de raiz média dos quadrados (RMS) que discutimos na segunda unidade.
• A raiz média dos quadrados dos erros (RMSE, do inglês root mean squared error) é defida por: \[\mbox{RMSE}=\sqrt{\frac{1}{n}SSR}=\sqrt{\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}^2}\]
• O RMSE representa a magnitude média dos erros de previsão da regressão.

• No R podemos encontrar o RMSE a partir dos resíduos.

epsilon.hat <- resid(fit)
sqrt(mean(epsilon.hat^2))

## [1] 0.2642361

• O resultado sugere que a medida de competência ajuda a prever o resultado da eleição, mas a previsão não é acurada. As previsões feitas com essa medida apresentam, em média, um erro de previsão de 26%.

• Os valores dos parâmetros obtidos com o método dos mínimos quadrados são dados por: \[\hat{\alpha}=\bar{Y}-\hat{\beta}\bar{X}\] \[\hat{\beta}=\frac{\sum_{i=1}^n(Y_i-\bar{Y})(X_i-\bar{X})}{\sum_{i=1}^n(X_i-\bar{X})^2}\]
• \(\bar{Y}=\frac{1}{n}\sum_{i=1}^nY_i\) e \(\bar{X}=\frac{1}{n}\sum_{i=1}^nX_i\) são as médias amostrais de \(Y\) e de \(X\).

• Note que a reta de regressão sempre passa por \((\bar{X},\bar{Y})\). Para ver que isso é verdade use a reta de regressão para calcular \(\hat{Y}\) associado a \(\bar{X}\) e substitua o valor de \(\hat{\alpha}\): \[\hat{Y}=\hat{\alpha}+\hat{\beta}\bar{X}=(\bar{Y}-\hat{\beta}\bar{X})+\hat{\beta}\bar{X}=\bar{Y}\]

• Uma propriedade importante do estimador de mínimos quadrados é que a média dos resíduos é zero, o que quer dizer que, em média, as previsões são acuradas. \[\mbox{média de }\hat{\varepsilon}=\frac{1}{n}\sum_{i=1}^n(Y_i-\hat{\alpha}-\hat{\beta}X_i)=\bar{Y}-\hat{\alpha}-\hat{\beta}\bar{X}=0\]
• O fato da média dos resíduos de uma regressão linear ser zero não significa necessariamente que o modelo estimado por mínimos quadrados é uma representação acurada do modelo gerador de dados.

mean(resid(fit))

## [1] -2.04732e-17

• O método dos mínimos quadrados é uma forma comum para estimar os parâmetros de um modelo de regressão linear. O método minimiza a soma dos quadrados dos resíduos, \(SSR=\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}^2=\frac{1}{n}\sum_{i=1}^n(Y_i-\hat{\alpha}-\hat{\beta}X_i)^2\). A média dos resíduos será sempre zero.

• É útil conhecer a relação entre a inclinação da reta de regressão e o coeficiente de correlação entre \(X\) e \(Y\). Para isso multiplique e divida \(\hat{\beta}\) pelo desvio padrão de Y, \(\sqrt{\frac{1}{n}\sum_{i=1}^n(Y_i - \bar{Y})^2}\): \[\hat{\beta}=\frac{1}{n}\sum_{i=1}^n\frac{(Y_i-\bar{Y})(X_i-\bar{X})}{\sqrt{\frac{1}{n}\sum_{i=1}^n(Y_i - \bar{Y})^2}\sqrt{\frac{1}{n}\sum_{i=1}^n(X_i - \bar{X})^2}}\\\times\frac{\sqrt{\frac{1}{n}\sum_{i=1}^n(Y_i - \bar{Y})^2}}{\sqrt{\frac{1}{n}\sum_{i=1}^n(X_i - \bar{X})^2}}\]

• Pelas definições de correlação e desvio padrão: \[\hat{\beta}=\mbox{correlação entre X e Y}\times\frac{\mbox{desvio padrão de Y}}{\mbox{desvio padrão de X}}\]

• A expressão acima tem duas implicações importantes:
  • A inclinação tem o mesmo sinal da correlação.
  • Um aumento de um desvio padrão em \(X\) é associado com um crescimento médio de \(\rho\) desvios padrão de \(Y\), onde \(\rho\) é a correlação entre \(X\) e \(Y\).

• No exemplo, a correlação entre a medida de competência, \(X\), e a diferença das margens, \(Y\), é de 0,46, o desvio padrão de \(X\) é 0,19 e o de \(Y\) é 0,29. Desta forma, um aumento de 0,19 na medida de competência está associado a um aumento médio de \(0,\!43 \times 0,\!29 = 0,\!1247\), ou 12,47%, na diferença das margens.

• Em 1886, Sir Francis Galton publicou um artigo chamado Regression towards mediocrity in herediatary stature, o artigo está entre os pioneiros no uso de análise de regressão. O trabalho estudava a relação entre a altura média dos país e dos filhos quando adultos.
• Galton foi o primeiro a observar um fenômeno conhecido como regressão à média (regression toward the mean). A conclusão de Galton foi de que quando a altura média dos parentes era menor que a média os filhos tendiam a ser mais altos do que a média.
• O fênomeno de regressão à média aparece em variáveis relacionadas à ciências sociais e economia.

• A regressão à média representa um fenômeno empírico onde uma observação com o valor do preditor muito distante da média da distribuição tende a ter um valor da variável de resultado mais próxima da média. Essa tendência pode ser explicada pelo acaso.

• Vamos avaliar a existência de regressão á média nas eleições americanas, mas, antes, é útil discutirmos sobre como juntar bases de dados no R.
• Fizemos isso usando cbind(), mas era preciso que as observações tivessem na mesma ordem, isso nem sempre é simples de conseguir. O tidyverse tem uma série de funções para juntar dois data.frames, \(df_1\) e \(df_2\):
  • left_join(): inclui todas as linhas do \(df_1\)
  • right_join(): inlui todas as linhas do \(df_2\)
  • full_join(): inclui as linhas comuns em \(df_1\) e \(df_2\)
  • inner_join(): inclui todas as linhas em \(df_1\) ou \(df_2\)

• Os dados em df1 e df2 mostram os cinco melhores classificados no Brasileirão de 2021 e 2022, respectivamente.

head(df1)

##               time pontos
## 1 Atlético Mineiro     84
## 2         Flamengo     71
## 3        Palmeiras     66
## 4        Fortaleza     58
## 5      Corinthians     57

head(df2)

##            time pontos
## 1     Palmeiras     81
## 2 Internacional     73
## 3    Fluminense     70
## 4   Corinthians     65
## 5      Flamengo     62

left_join(df1,df2, by=c("time"="time"), suffix=c(".21",".22"))

##               time pontos.21 pontos.22
## 1 Atlético Mineiro        84        NA
## 2         Flamengo        71        62
## 3        Palmeiras        66        81
## 4        Fortaleza        58        NA
## 5      Corinthians        57        65

right_join(df1,df2, by=c("time"="time"), suffix=c(".21",".22"))

##            time pontos.21 pontos.22
## 1      Flamengo        71        62
## 2     Palmeiras        66        81
## 3   Corinthians        57        65
## 4 Internacional        NA        73
## 5    Fluminense        NA        70

inner_join(df1,df2, by=c("time"="time"), suffix=c(".21",".22"))

##          time pontos.21 pontos.22
## 1    Flamengo        71        62
## 2   Palmeiras        66        81
## 3 Corinthians        57        65

full_join(df1,df2, by=c("time"="time"), suffix=c(".21",".22"))

##               time pontos.21 pontos.22
## 1 Atlético Mineiro        84        NA
## 2         Flamengo        71        62
## 3        Palmeiras        66        81
## 4        Fortaleza        58        NA
## 5      Corinthians        57        65
## 6    Internacional        NA        73
## 7       Fluminense        NA        70

• Começamos nossa análise comparando a fração de votos obtida por Obama em 2008 e 2012. Os dados referentes às eleições de 2012 estão no data.frame pres12.
• O primeiro passo é juntar os dados em pres08 e pres12. Usaremos a função full_join().

• Variáveis em pres12:
  • state: sigla do estado
  • Obama: percentual dos votos recebido por Obama
  • Rommey: percentual dos votos recebido por Rommey
  • EV: número de votos do estado no colégio eleitoral

pres12 <- read.csv("pres12.csv")
head(pres12)

##   state Obama Romney EV
## 1    AL    38     61  9
## 2    AK    41     55  3
## 3    AZ    45     54 11
## 4    AR    37     61  6
## 5    CA    60     37 55
## 6    CO    51     46  9

pres <- full_join(x=pres08, y=pres12, by="state")

str(pres)

## 'data.frame':    51 obs. of  11 variables:
##  $ state.name: chr  "Alabama" "Alaska" "Arizona" "Arkansas" ...
##  $ state     : chr  "AL" "AK" "AZ" "AR" ...
##  $ Obama.x   : int  39 38 45 39 61 54 61 92 62 51 ...
##  $ McCain    : int  60 59 54 59 37 45 38 7 37 48 ...
##  $ EV.x      : int  9 3 10 6 55 9 7 3 3 27 ...
##  $ margin    : int  -21 -21 -9 -20 24 9 23 85 25 3 ...
##  $ poll.pred : num  -25 -19 -2.5 -7 24 7 25 69 30 2 ...
##  $ correct   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Obama.y   : int  38 41 45 37 60 51 58 91 59 50 ...
##  $ Romney    : int  61 55 54 61 37 46 41 7 40 49 ...
##  $ EV.y      : int  9 3 11 6 55 9 7 3 3 29 ...

• Como exercício vamos melhorar os nomes das variáveis comuns aos dois data.frames e mudar o nome da coluna com as siglas dos estados (coluna usada para juntar as bases de dados) em pres12.

pres12 <- rename(pres12, state.abbrev=state)
head(pres12)

##   state.abbrev Obama Romney EV
## 1           AL    38     61  9
## 2           AK    41     55  3
## 3           AZ    45     54 11
## 4           AR    37     61  6
## 5           CA    60     37 55
## 6           CO    51     46  9

pres <- full_join(pres08,pres12, by = c("state" = "state.abbrev"),
                  suffix = c("_08", "_12"))

• Por conta da polarização crescente na política americana, antes de fazer as comparações entre a votação em 2008 e 2012 os percentuais de votos em Obama serão padronizados, ou seja, vamos usar o valor z nas variáveis com o percentual de votos em Obama em 2008 e 2012.
• A função scale() faz esse trabalho.

pres <- pres %>%
  mutate(Obama2008.z = as.numeric(scale(Obama_08)),
         Obama2012.z = as.numeric(scale(Obama_12)))

• A regressão será entre os votos padronizados para Obama em 2012 e 2008. Esperamos uma forte relação positiva entre os dois, os votos para Obama em 2012 tendem a ser altos nos estados onde Obama foi bem votado em 2008.
• Como a regressão é entre duas variáveis padronizadas o intercepto será zero. Isso é verdade poruqe, devido a padronização \(\bar{Y}\) e \(\bar{X}\) são iguais a zero e \(\hat{\alpha} = \bar{Y} - \hat{\beta}\bar{X}\).
• É possível fazer a regressão sem a constante adicionando \(-1\) na fórmula.

fit1 <- lm(Obama2012.z ~ Obama2008.z, data = pres)
fit1

## 
## Call:
## lm(formula = Obama2012.z ~ Obama2008.z, data = pres)
## 
## Coefficients:
## (Intercept)  Obama2008.z  
##  -5.076e-17    9.834e-01

fit1 <- lm(Obama2012.z ~ -1 + Obama2008.z, data = pres)
fit1

## 
## Call:
## lm(formula = Obama2012.z ~ -1 + Obama2008.z, data = pres)
## 
## Coefficients:
## Obama2008.z  
##      0.9834

• Vamos “dar uma olhada” na regressão.

pres %>%
  ggplot(aes(Obama2008.z,Obama2012.z)) +
  geom_smooth(method="lm", se=FALSE, linewidth=0.5) +
  geom_point(shape=1) +
  coord_fixed() +
  scale_x_continuous(limits = c(-4,4)) +
  scale_y_continuous(limits = c(-4,4)) +
  labs(x = "Obama's standardized votes share in 2008",
       y = "Obama's standardized votes share in 2012") +
  theme_classic()

• Para avaliar se há regressão à média vamos calcular a proporção de estados onde Obama recebeu uma fração de votos maior em 2012 do que em 2008. Primeiro vamos fazer a conta para os estados no quartil inferior dos votos padronizados em Obama nas eleições de 2008, depois para o quartil superior.
• A regressão à média será observada se a proporção for maior no quartil inferior do que no quartil superior.

pres %>%
  filter(Obama2008.z < quantile(Obama2008.z, 0.25)) %>%
  summarize(improve = mean(Obama2012.z > Obama2008.z))

##     improve
## 1 0.5833333

pres %>%
  filter(Obama2008.z > quantile(Obama2008.z, 0.75)) %>%
  summarize(improve = mean(Obama2012.z > Obama2008.z))

##     improve
## 1 0.4615385

• Obama melhorou o desempenho em 58% dos estados do quartil inferior em 2008 (grupo onde teve pior desempenho), no quartil superior (onde teve melhor desempenho), Obama melhorou em apenas 46% dos estados.

• Quão bem um modelo se ajusta aos dados? Existem medidas para responder essa pergunta. A mais conhecida é o coeficiente de determinação, ou \(R^2\).
• O \(R^2\) representa a proporção da variação total de \(Y\) que é explicada pelo modelo.

• Para definir o \(R^2\) primeiro vamos definir a soma total dos quadrados (TSS, do inglês: total sum of squares): \[\mbox{TSS}=\sum_{i=1}^n(Y_i - \bar{Y})^2\]
• O \(R^2\) é a proporção do TSS explicada pelo preditor \(X\): \[R^2=1-\frac{\mbox{SSR}}{\mbox{TSS}}\]

• Se a soma dos quadrados dos resíduos, SSR, que é a variação não explicada pelo modelo, for próxima de TSS, então \(\frac{\mbox{SSR}}{\mbox{TSS}}\) será próxima de 1 e o \(R^2\) próximo de zero. O modelo se ajusta pouco aos dados.
• Por outro lado, se SSR for pequena em relação a TSS, então \(\frac{\mbox{SSR}}{\mbox{TSS}}\) será próxima de zero e o \(R^2\) será próximo de um. O modelo se ajusta bem aos dados.

• O coeficiente de determinação, também chamado \(R^2\), é uma medida de ajuste do modelo aos dados e representa a proporção da variação de \(Y\) que é explicada por \(X\). É definido como um menos a razão entre a soma dos quadrados dos resíduos e a soma total dos quadrados.

• Como exemplo considere um modelo para prever o resultado das eleições presidenciais na Flórida em 2000 usando dados a nível de condados referentes as eleições de 1996.
• Os dados estão no data.frame florida.

• Variáveis:
  • county: nome do condado
  • Clinton96: votos para Clinton em 1996
  • Dole96: votos para Dole em 1996
  • Perot96: votos para Perot em 1996
  • Bush00: votos para Bush em 2000
  • Gore00: votos para Gore em 2000
  • Buchanan00: votos para Buchanan em 2000

florida <- read.csv("florida.csv")

• Vamos focar nos candidatos do Partido Reformista, Reform Party, um partido fundado por Ross Perot em 1995.

fit2 <- lm(Buchanan00 ~ Perot96, data=florida)
fit2

## 
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida)
## 
## Coefficients:
## (Intercept)      Perot96  
##     1.34575      0.03592

• Vamos calcular o TSS, o SSR e o \(R^2\):

TSS2 <- florida %>%
  mutate(diff_sq = (Buchanan00 - mean(florida$Buchanan00))^2) %>%
  summarize(TSS=sum(diff_sq))

SSR2 <- sum(resid(fit2)^2)

1 - SSR2/TSS2

##         TSS
## 1 0.5130333

• O resultado diz que os votos de Perot em 1996 explicam 51% da variação dos votos de Buchanan em 2000. Que tal criar uma função para calcular o \(R^2\)?

R2 <- function(fit){
  resid <- resid(fit)
  y <- fitted(fit) + resid
  TSS <- sum((y-mean(y))^2)
  SSR <- sum(resid^2)
  R2 <- 1 - SSR/TSS
  return(R2)
}

R2(fit2)

## [1] 0.5130333

• A função summary() quando aplicada a um objeto do tipo lm retorna informações relativas a regressão, entre essas informações está o \(R^2\).

summary(fit2)

## 
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -612.74  -65.96    1.94   32.88 2301.66 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.34575   49.75931   0.027    0.979    
## Perot96      0.03592    0.00434   8.275 9.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 316.4 on 65 degrees of freedom
## Multiple R-squared:  0.513,  Adjusted R-squared:  0.5055 
## F-statistic: 68.48 on 1 and 65 DF,  p-value: 9.474e-12

• Caso o interesse seja apenas o \(R^2\) podemos usar a função summary() e buscar o \(R^2\):

fit2summary <- summary(fit2)
fit2summary$r.squared

## [1] 0.5130333

• Também podemos conseguir o \(R^2\) usando a função glance() do pacote broom.

glance(fit2)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.513         0.506  316.      68.5 9.47e-12     1  -480.  966.  972.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

• Vimos que os votos de Obama em 2008 tem uma correlação forte com os votos de Obama em 2012, isso poderia nos fazer esperar um \(R^2\) maior do que o que encontramos na regressão para os candidatos do Partido Reformista na Flórida.
• O \(R^2\) da regressão entre votos em Obama em 2008 e 2012 foi de 97%!

R2(fit1)

## [1] 0.9671579

• Olhar os dados é sempre um bom caminho para entender “resultados estranhos”. Vamos fazer um gráfico de dispersão com os resíduos e o s valores ajustados para variável dependente.

augment_fit2 <- augment(fit2)

augment_fit2 %>%
  ggplot(aes(.fitted, .resid)) +
  geom_point(shape=1) +
  geom_hline(yintercept = 0) +
  labs(x="Fitted values", y="Residuals") +
  theme_classic()

• Um dos resíduos é extremamente alto, nesses casos falamos de um outlier, nesse condado Buchanan teve mais de 2000 votos além do esperado pelo modelo. Que condado foi esse?
• As funções add_predictions() e add_residuals() do pacote modelr ajudam a identificar o condado.

library(modelr)

florida_fit2 <- florida %>%
  add_predictions(fit2) %>%
  add_residuals(fit2)

florida_fit2 %>%
  filter(resid == max(resid)) %>%
  select(county) %>%
  pull()

## [1] "PalmBeach"

• Uma possível explicação para o comportamento estranho dos votos em Palm Beach é adoção da cédula do tipo borboleta por este condado. Nesse tipo de cédula os eleitores devem perfurar no ponto correspondente ao candidato que desejam votar, mas o formato da cédula pode induzir eleitores ao erro.

• Bush ganhou na Flórida por uma margem de 537 votos, os votos da Flórida no colégio eleioral decidiram as eleições de 2000.
• Muita gente acredita que anomalias em Palm Beach foram responsáveis pela derrota de Gore.
• Para terminar o exemplo, façamos o exercício navamente exluindo o condado de Palm Beach.

florida.pb <- filter(florida, county != "PalmBeach")

fit3 <- lm(Buchanan00 ~ Perot96, data = florida.pb)
fit3

## 
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida.pb)
## 
## Coefficients:
## (Intercept)      Perot96  
##    45.84193      0.02435

R2(fit3)

## [1] 0.8511675

• O \(R^2\) foi de 0,513 para 0,851.
• Na sequência vamos fazer o gráfico dos resíduos e um gráfico com a reta de regressão com Palm Beach e a reta de regressão sem Palm Beach.

florida.pb %>%
  add_residuals(fit3) %>%
  add_predictions(fit3) %>%
  ggplot(aes(pred, resid)) +
  geom_hline(yintercept = 0) +
  geom_point(shape=1) +
  ylim(-750, 2500) +
  xlim(0,1500) +
  labs(title = "Residual plot without Plam Beach",
       x= "Fitted values",
       y="Residuals") +
  theme_classic()

ggplot() +
  geom_point(data = florida, aes(Perot96, Buchanan00), shape=1) +
  geom_smooth(data = florida, aes(Perot96, Buchanan00), method="lm", se=FALSE) +
  geom_smooth(data = florida.pb, aes(Perot96, Buchanan00), method="lm", 
              se=FALSE, linetype="dashed") +
  annotate("text", x=30000, y=3200, label = "Palm Beach") +
  annotate("text", x=30000, y=300, label = "Regression line without\n Palm Beach") +
  annotate("text", x=25000, y=1400, label = "Regression line with\n Palm Beach") +
  theme_classic()

• Os modelos que estudamos nesta seção tem por base previsões dentro da amostra (in-sample predictions) em vez de previsões fora da amostra (out-of-sample predictions). Isso quer dizer que as estatística do modelo, como o \(R^2\), descrevem quão bem o modelo se ajusta a amostra usada para estimar o modelo. Se o modelo for desenhado para ajustar muito bem em uma determinada amostra, quando isso ocorre falamos de overfitting, o modelo pode fazer um ajuste ruim para outras amostras.
• Quando desejamos um modelo que possa ser generalizado para outros dados, temos de tomar cuidado para evitar overfitting do modelo para uma amostra específica.

• Análise de regressão é uma ferramenta básica para previsão em economia, boa parte dos cursos de econometria é dedicada a formas de estimar os parâmetros da reta de regressão que melhor se adaptem as características dos dados que usamos em economia.
• Porém, é preciso ter claro que apenas em condições bem específicas regressão pode ser usada para estabelecer inferência causal. A regressão quantifica relação entre variáveis e, como sabemos, relação não implica causalidade.
• Para falar de causalidade é preciso ter o contrafactual, em alguns casos é possível prever o contrafactual usando regressão. Fora desses casos, falar de causalidade com resultados de regressão é um erro que pode ser perigoso.

• Voltemos aos experimentos aleatórios, mas agora sabendo fazer regressões.
• O exemplo é um estudo que buscou avaliar o efeito causal nas políticas públicas de mulheres no governo. Especificamente, o estudo quer saber se mulheres no governo promovem políticas diferentes das promovidas por homens.
• O exemplo segue o artigo: Raghabendra Chattopadhyay e Esther Duflo. Women as policy makers: evidence from a ramdomized experimento in India. Econometrica, v.72, n.5, 2004.

• Uma comparação entre distritos que elegem homens e distritos que elegem mulheres não é adequada para responder a questão do exemplo. É possível que existam outros fatores (confundidores) que associem o tipo de política a eleição de mulheres. Por exemplo, é possível que distritos mais liberais, onde os eleitos tendem a aplicar políticas mais liberais, tendam a eleger mais mulheres.
• Para driblar esse potencial confundidor os autores do estudo se aproveitaram de um experimento aleatório que aconteceu na Índia em meados da década de 1990: um terço das chefias de conselho da vilas foram aleatoriamente reservado para mulheres.

• A política foi implementada em um nível de governo chamado Gram Panchayat, ou GP. cada GP contém várias vilas.
• Os dados em women correspondem ao estado de West Bengal. Duas vilas de cada GP foram selecionada aleatoriamente para coleta detalhada de dados.
• Cada observação na base de dados corresponde a uma vila, são duas vilas de cada GP.

• Variáveis:
  • GP: indetificador para Gram Pachayat (GP)
  • village: identificador para vila
  • reserved: variável binária indicando se o GP tem líder mulher ou não
  • irrigation: valiável medindo o número de instalações para irrigação (novas ou reformadas) na vila desde que a política de reservas foi iniciada.
  • water: variável medindo o número de instalações para água para beber (novas ou reformadas) na vila desde que a política de reservas foi iniciada.

women <- read.csv("women.csv")
head(women)

##   GP village reserved female irrigation water
## 1  1       2        1      1          0    10
## 2  1       1        1      1          5     0
## 3  2       2        1      1          2     2
## 4  2       1        1      1          4    31
## 5  3       2        0      0          0     0
## 6  3       1        0      0          0     0

• A primeira tarefa é checar se a lei foi seguida, para isso vamos calcular a proporção de mulheres eleitas para os cargos reservados e não reservados. Como cada GP tem o mesmo número de vilas podemos simplesmente calcular as médias das vilas sem criar uma base de dados a nível de GP.
• Se a lei foi cumprida as proporção nos cargos reservados deve ser um.

women %>%
  group_by(reserved) %>%
  summarize(prop_female = mean(female))

## # A tibble: 2 × 2
##   reserved prop_female
##      <int>       <dbl>
## 1        0      0.0748
## 2        1      1

• Os dados mostram que a lei foi seguida. Em todos os GPs que deviam reservar vagas para mulheres há pelo menos uma mulher eleita. Nos que não tinham de reservar vagas para mulheres apenas 7,5% tem pelo menos uma mulher eleita, ou seja, em 92,5% dos GPs que não tinham reservas para mulheres não há mulher eleita.
• Como a dsitribuição dos GPs com reserva para mulheres foi aleatória os GPs sem essa reserva formam um contrafactual.

• A hipótese a ser testada é se políticas mulheres estão mais inclinadas a atender demandas das mulheres eleitoras.
• Pesquisas mostram que mulheres reclamam mais frequentemente sobre a qualidade da água para beber e homens reclamam com mais frequência de falta de irrigação.
• Vamos estimar o efeito causal da reserva de vagas para mulheres no número de estrutras de irrigação e água para beber, nos dois casos serão consideradas estruturas novas ou reformadas. A estimação será feita por meio de diferenças em médias (vimos esse estimador na segunda unidade).

women %>%
  group_by(reserved) %>%
  summarize(irrigation = mean(irrigation),
            water = mean(water)) %>%
  pivot_longer(names_to = "variables", -reserved) %>%
  pivot_wider(names_from = reserved) %>%
  rename("not_reserved" = `0`, "reserved" = `1`) %>%
  mutate(diff = reserved - not_reserved)

## # A tibble: 2 × 4
##   variables  not_reserved reserved   diff
##   <chr>             <dbl>    <dbl>  <dbl>
## 1 irrigation         3.39     3.02 -0.369
## 2 water             14.7     24.0   9.25

• Os resultados mostram um efeito positivo das reservas em estruturas de água para beber, o efeito nas estruturas de irrigação foi negativo, porém bem pequeno.
• Os resultado indicam que a reserva para mulheres de fato aumentou o atendimento de demandas feitas por mulheres.

• Podemos usar de regressões para obter esse mesmo resultado.
• Uma regressão entre a variável de resultado em uma variável de tratamento gera uma reta onde a inclinação é igual a diferença na variável de resultado entre os dois grupos. O intercepto corresponde a média da variável de resultado no grupo de controle.

• De forma mais geral, quando \(X\) é uma variável binária os coeficientes estimados em uma regressão são: \[\hat{\alpha}=\underbrace{\frac{1}{n_0}\sum_{i=1}^n(1-X_i)Y_i}_{\mbox{média do grupo de controle}}\] \[\hat{\beta}=\underbrace{\frac{1}{n_1}\sum_{i=1}^nX_iY_i}_{\mbox{média do grupo de tratamento}} - \underbrace{\frac{1}{n_0}\sum_{i=1}^n(1-X_i)Y_i}_{\mbox{média do grupo de controle}}\]
• Onde \(n_0\) e \(n_1\) representam o número de observações nos grupos de controle e tratamento.

• Ver para crer:

lm(water ~ reserved, data = women)

## 
## Call:
## lm(formula = water ~ reserved, data = women)
## 
## Coefficients:
## (Intercept)     reserved  
##      14.738        9.252

lm(irrigation ~ reserved, data = women)

## 
## Call:
## lm(formula = irrigation ~ reserved, data = women)
## 
## Coefficients:
## (Intercept)     reserved  
##      3.3879      -0.3693

• Para “enxergar” esse resultado lembre que a variável \(X\) só tem valores zero e um, zero para o grupo de controle e um para o grupo tratado.
• O valor estimado para \(Y\) quando \(X=0\) é: \[\hat{Y(0)}=\hat{\alpha} + \hat{\beta}\times0= \hat{\alpha}\]
• O valor estimado para diferença de \(Y\) quando \(X=1\) e quando \(X=0\) é: \[\hat{Y(1)}-\hat{Y(0)}=\hat{\alpha} + \hat{\beta}\times1 - (\hat{\alpha} + \hat{\beta}\times0)=\hat{\alpha} + \hat{\beta}-\hat{\alpha}=\hat{\beta}\]

• Quando aplicado em dados experimentais, o valor estimado para inclinação de um modelo de regressão linear pode ser interpretado como uma estimativa do efeito médio do tratamento é numericamente equivalente ao estimador de diferenças em médias. O valor estimado para o intercepto é igual a média da variável de resultado para o grupo de controle. A distribuição do tratamento ser aleatória que permite a interpretação desses estimadores como um efeito causal.

• Até agora discutimos regressão com apenas uma variável explicativa, o chamado modelo de regressão linear simples, mas esse não precisa ser o caso. MUitas vezes o modelo tem mais de uma variável explicativa, quando isso ocorre falamos de regressão linear múltipla que tem a seguinte forma: \[Y=\alpha+\beta_1 X_1+\beta_2 X_2 + \cdots + \beta_p X_p + \varepsilon\]
• No modelo \(\alpha\) é o intercepto, \(\beta_j\) é coeficiente da variável \(X_j\) e \(\varepsilon\) é o termo de erro. O modelo tem \(p\) variáveis explicativas.

• No modelo de regressão linear múltipla cada \(\beta_j\) representa a variação em \(Y\) associada ao aumento de uma unidade em \(X_j\) quando todas as outras variáveis explicativas estão constantes, é o que chamamos de ceteris paribus.
• Desta forma, o modelo de regressão múltipla permite avaliar o impacto em \(Y\) de cada variável explicativa.

• O método dos mínimos quadrados pode ser utilizado para obter os estimadores do intercepto e de cada coeficiente do modelo de regressão múltipla.
• A soma dos quadrados dos resíduos, SSR, é dada por: \[\mbox{SSR}=\sum_{i=1}^n\hat{\varepsilon}=\sum_{i=1}^n(Y_i - \hat{\alpha} - \hat{\beta}_1 X_{i1} - \hat{\beta}_2 X_{i2} - \cdots - \hat{\beta}_p X_{ip})^2\]

• Como exemplo de regressão múltipla considere o estudo de pressão social e comparecimento nas eleições que vimos na Unidade 2.
• Para fazer regressões múltiplas usamos a função lm(), a fórmula no chamado da função deve incluir todas as variáveis explicativas, por exemplo: lm(y ~ x1 + x2 + x3).
• No modelo do exemplo a variável explicativa, messages, é uma variável da classe character com os vários grupos.Nesse casos a função lm() automaticamente cria um conjunto de variáveis binárias com um nas observações correspondentes ao grupo e zero nas outras observações.

social <- read.csv("social.csv")
class(social$messages)

## [1] "character"

head(social$messages)

## [1] "Civic Duty" "Civic Duty" "Hawthorne"  "Hawthorne"  "Hawthorne" 
## [6] "Control"

unique(social$messages)

## [1] "Civic Duty" "Hawthorne"  "Control"    "Neighbors"

fit <- lm(primary2006 ~ messages, data=social)
fit

## 
## Call:
## lm(formula = primary2006 ~ messages, data = social)
## 
## Coefficients:
##       (Intercept)    messagesControl  messagesHawthorne  messagesNeighbors  
##          0.314538          -0.017899           0.007837           0.063411

• Se preferir, você pode criar as variáveis indicadoras de cada grupo e fazer a regressão.

social <- social %>%
  mutate(Control = if_else(messages == "Control", 1, 0),
         Hawthorne = if_else(messages == "Hawthorne", 1, 0),
         Neighbors = if_else(messages == "Neighbors", 1, 0))

social %>%
  select(Control, Hawthorne, Neighbors) %>%
  head()

##   Control Hawthorne Neighbors
## 1       0         0         0
## 2       0         0         0
## 3       0         1         0
## 4       0         1         0
## 5       0         1         0
## 6       1         0         0

lm(primary2006 ~ Control + Hawthorne + Neighbors, data = social)

## 
## Call:
## lm(formula = primary2006 ~ Control + Hawthorne + Neighbors, data = social)
## 
## Coefficients:
## (Intercept)      Control    Hawthorne    Neighbors  
##    0.314538    -0.017899     0.007837     0.063411

• Repare que o grupo Civic Duty não apareceu na regressão, isso acontece porque o grupo foi escolhido como grupo de referência.
• Dessa forma a proproção de eleitores que compareceram no grupo Civic Duty é dado pela constante (pois as variáveis dos outros grupos são iguais a zero) e cada coeficiente \(\beta\) representa a diferença da proporção dos eleitores que compareceram no grupo e essa mesma proporção no grupo Civic Duty.

• Por exemplo, no grupo Control o valor previsto para \(Y\) será: \(\hat{\alpha} + \hat{\beta}_1 = 0,\!315 + (-0,\!018) = 0,\!297 = 29,7\%\)
• Do mesmo modo, para o grupo Neighbor o valor previsto para \(Y\) será: \(\hat{\alpha} + \hat{\beta}_3 = 0,\!315 + 0,\!063 = 0,\!378 = 37,8\%\)

• Podemos obter todos os valores previstos usando a função add_predictions() do pacote modelr. Para isso vamos usar a função data_grid(), também do modelr, para criar um data.frame com cada valor único da variável messages. A função data_grid() cria um data.frame com combinações únicas das colunas especificadas.

data_grid(social, messages)

## # A tibble: 4 × 1
##   messages  
##   <chr>     
## 1 Civic Duty
## 2 Control   
## 3 Hawthorne 
## 4 Neighbors

unique_messages <- data_grid(social, messages) %>%
  add_predictions(fit)

unique_messages

## # A tibble: 4 × 2
##   messages    pred
##   <chr>      <dbl>
## 1 Civic Duty 0.315
## 2 Control    0.297
## 3 Hawthorne  0.322
## 4 Neighbors  0.378

• Assim como na regressão simples, os valores previstos na regressão múltipla correspondem as médias de cada grupo.

social %>%
  group_by(messages) %>%
  summarize(mean(primary2006))

## # A tibble: 4 × 2
##   messages   `mean(primary2006)`
##   <chr>                    <dbl>
## 1 Civic Duty               0.315
## 2 Control                  0.297
## 3 Hawthorne                0.322
## 4 Neighbors                0.378

• Caso você queira obter os valores previstos direto da regressão, em vez de obter o valor previsto para o grupo de referência e as diferenças em relação a esse grupo, basta fazer a regressão sem a constante.

fit.noint <- lm(primary2006 ~ -1 + messages, data=social)
fit.noint

## 
## Call:
## lm(formula = primary2006 ~ -1 + messages, data = social)
## 
## Coefficients:
## messagesCivic Duty     messagesControl   messagesHawthorne   messagesNeighbors  
##             0.3145              0.2966              0.3224              0.3779

• Usando fatores é relativamente simples mudar o grupo de referência, basta usar a função relevel() para mudar o nível de referência do fator.

social.f <- social %>%
  mutate(messages.f = as.factor(messages),
         mes.C = relevel(messages.f, ref = "Control"),
         mes.H = relevel(messages.f, ref = "Hawthorne"),
         mes.N = relevel(messages.f, ref = "Neighbors"))

lm(primary2006 ~ mes.C, data = social.f)

## 
## Call:
## lm(formula = primary2006 ~ mes.C, data = social.f)
## 
## Coefficients:
##     (Intercept)  mes.CCivic Duty   mes.CHawthorne   mes.CNeighbors  
##         0.29664          0.01790          0.02574          0.08131

lm(primary2006 ~ mes.H, data = social.f)

## 
## Call:
## lm(formula = primary2006 ~ mes.H, data = social.f)
## 
## Coefficients:
##     (Intercept)  mes.HCivic Duty     mes.HControl   mes.HNeighbors  
##        0.322375        -0.007837        -0.025736         0.055574

lm(primary2006 ~ mes.N, data = social.f)

## 
## Call:
## lm(formula = primary2006 ~ mes.N, data = social.f)
## 
## Coefficients:
##     (Intercept)  mes.NCivic Duty     mes.NControl   mes.NHawthorne  
##         0.37795         -0.06341         -0.08131         -0.05557

• Se o interesse é apenas nos valores das médias, você pode usar a função group_by() e calcular a diferença em relação ao grupo de interesse.

social %>%
  group_by(messages) %>%
  summarize(primary2006 = mean(primary2006)) %>%
  mutate(Control = primary2006[messages == "Control"],
         diff = primary2006 - Control)

## # A tibble: 4 × 4
##   messages   primary2006 Control   diff
##   <chr>            <dbl>   <dbl>  <dbl>
## 1 Civic Duty       0.315   0.297 0.0179
## 2 Control          0.297   0.297 0     
## 3 Hawthorne        0.322   0.297 0.0257
## 4 Neighbors        0.378   0.297 0.0813

• Por fim, podemos calcular o \(R^2\) para modelos de regressão múltipla, porém é comum usar o \(R^2\)-*, que ajusta o \(R^2\) pelos graus de liberdade, uma forma de corrigir pelo número de variáveis explicativas.
• Os grau de liberdade é dado pela diferença entre o número de observações, \(n\), e o número de parâmetros estimados, \(p+1\). Desta forma o grau de liberdade é igual a \(n-(p+1)\).

• O \(R^2\) ajustado é dado por: \[R^2-\mbox{ajustado}=1-\frac{\mbox{SSR}/(n-p-1)}{TSS/(n-1)}\]

• Assim como fizemos para o \(R^2\), podemos criar uma função para calcular o \(R^2\)-ajustado.

adjR2 <- function(fit) {
  resid <- resid(fit)
  y <- fitted(fit) + resid
  n <- length(y)
  TSS.adj <- sum((y - mean(y))^2)/(n - 1)
  SSR.adj <- sum(resid^2)/(n - length(coef(fit)))
  R2.adj <- 1 - SSR.adj/TSS.adj
  return(R2.adj)
}

adjR2(fit)

## [1] 0.003272788

R2(fit)

## [1] 0.003282564

summary(fit)$adj.r.squared

## [1] 0.003272788

• O modelo de regressão linear com múltiplas variáveis explicativas é definido por: \[Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \varepsilon\] Onde o coeficiente \(\beta_j\) representa o aumento médio da variável de resultado associado a aumento de uma unidade em \(X_j\) mantendo as demais variáveis explicativas constantes. Os coeficientes são estimados pelo método de minimizar a soma dos quadrados dos resíduos. Quando do cálculo do \(R^2\) é comum ajustar pelos graus de liberdade.

• Quando aplicados em experimentos aleatórios, o modelo de regressão linear ajuda a explorar os efeitos do tratamento em diferentes grupos, os chamados efeitos heterogêneos do tratamento.
• O tratamento pode não ter o mesmo efeito em diferentes unidades tratadas, desta forma identificar características associadas a direção e magnitude do tratamento é essencial para determinar quem deve receber o tratamento.
• Quais perfis de famílias são mais beneficiadas por programas de transferência de renda? Quais tipos de empresas são mais afetadas por uma política de redução de barreiras à entrada?

• Para dar exemplo desse uso do modelo de regressão vamos retornar à questão da pressão social sobre os eleitores. Queremos agora saber se a pressão social tem mais (ou menos) efeito em eleitores que votam com frequência.
• Em termos práticos, vamos comparar o efeito causal da mensagem do tipo Neighbor nos eleitores que votaram e dos que não votaram em 2004.
• O procedimento será dividir os dados e calcular o efeito médio do tratamento (ate do inglês: average treatment effect) no grupo dos que votaram em 2004 e no grupo dos que não votaram em 2004.

ate <- social %>%
  group_by(primary2004, messages) %>%
  summarize(primary2006 = mean(primary2006)) %>%
  pivot_wider(names_from = messages, values_from = primary2006) %>%
  mutate(ate_Neighbors = Neighbors - Control) %>%
  select(primary2004, Neighbors, Control, ate_Neighbors) %>%
  ungroup()

ate

## # A tibble: 2 × 4
##   primary2004 Neighbors Control ate_Neighbors
##         <int>     <dbl>   <dbl>         <dbl>
## 1           0     0.306   0.237        0.0693
## 2           1     0.482   0.386        0.0965

ate.voter <- filter(ate, primary2004 == 1) %>%
  select(ate_Neighbors) %>% pull()

ate.nonvoter <- filter(ate, primary2004 == 0) %>%
  select(ate_Neighbors) %>% pull()

ate.voter - ate.nonvoter

## [1] 0.02722908

• O exercício mostra que o efeito das cartas sobre o comparecimento dos que votaram em 2004 foi de 9,65% e nos que não votaram em 2003 foi de 6,93%, a diferença foi de aproximadamente 2,7%.
• Portanto, as mensagens do tipo Neighbor afetaram bem mais os que participaram das eleições de 2004 do que os que não participaram.
• Podemos chegar a mesma conclusão usando um modelo de regressão linear com um efeito de interação entre a variável de tratamento Neighbors e variável explicativa primary2004.

• O modelo com o efeito de interação tem a forma: \[Y=\alpha+\beta_1\mbox{primary2004}+\beta_2\mbox{Neighbors}+ \\ +\beta_3(\mbox{primary2004}\times\mbox{Neighbors}) + \varepsilon\]
• Note que o produto entre primary2004 e Neighbors será igual a 1 se o indivíduo votou nas primárioas de 2004 e recebeu carta do tipo Neighbors, em todos os outros casos o produto será zero.

• Desta forma, entre os que compareceram às primárias de 2004 (primary2004=1) o efeito médio da mensagem do tipo Neighbors será de \(\beta_2 + \beta_3\). O efeito sobre os que não compareceram (primary2004=0) será igual a \(\beta_2\).
• Por fim, o coeficiente \(\beta_3\) mede a diferença do efeito entre os dois grupos.

• De forma mais geral, considere o exemplo de modelo de regressão linear com termo de interação: \[Y=\alpha+\beta_1 X_1+\beta_2 X2 + \beta_3 X_1 X_2 + \varepsilon\]
• O coeficiente \(\beta_3\) representa como o efeito de \(X_1\) depende de \(X_2\).

• Para ver que isso é verdade fixe \(X_2 = x_2\) e calcule o valor previsto de \(Y\) quando \(X_1 = x_1\). Esse valor será dado por: \(\hat{\alpha} + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \hat{\beta}_3 x_1 x_2\).
• Agora compare com o valor previsto quando \(X_1 = x_1 + 1\): \(\hat{\alpha} + \hat{\beta}_1 (x_1+1) + \hat{\beta}_2 x_2 + \hat{\beta}_3 (x_1+1) x_2\).
• Diminuindo o segundo valor do primeiro valor temos: \(\hat{\beta}_1+\hat{\beta}_3 X_2\).
• A expressão acima é uma equação linear. O intercepto, \(\beta_1\), representa o aumento médio na variável de resultado decorrente de um aumento de \(X_1\) quando \(X_2 = 0\). A inclinação, \(\beta_3\), representa o aumento médio em \(X_1\) quando \(X_2\) aumenta em uma unidade.

• A equação abaixo é um exemplo de modelo de regressão linear com termo de interação: \[Y=\alpha+\beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1X_2 + \varepsilon\] O modelo supõe que que o efeito de \(X_1\) depende linearmente de \(X_2\). Desta forma, quando \(X_2\) aumenta em uma unidade, a mudança em \(Y\) associada ao aumento de uma unidade de \(X_1\) é dada por \(\beta_3\).

• No R o termo de interação é representado por :, entra na fórmula como \(x_1:x_2\).

social.neighbor <- social %>%
  filter(messages == "Control" | messages == "Neighbors")

fit.int <- lm(primary2006 ~ primary2004 + messages + primary2004:messages,
              data = social.neighbor)

fit.int

## 
## Call:
## lm(formula = primary2006 ~ primary2004 + messages + primary2004:messages, 
##     data = social.neighbor)
## 
## Coefficients:
##                   (Intercept)                    primary2004  
##                       0.23711                        0.14870  
##             messagesNeighbors  primary2004:messagesNeighbors  
##                       0.06930                        0.02723

• Forma alternativa para regressão com termo de interação:

lm(primary2006 ~ primary2004*messages, data = social.neighbor)

## 
## Call:
## lm(formula = primary2006 ~ primary2004 * messages, data = social.neighbor)
## 
## Coefficients:
##                   (Intercept)                    primary2004  
##                       0.23711                        0.14870  
##             messagesNeighbors  primary2004:messagesNeighbors  
##                       0.06930                        0.02723

• Para interpretar o resultado da regressão lebre que o efeito médio do tratamento de receber uma carta do tipo Neighbor é a diferença no \(Y\) estimado no grupo de tratamento e controle.
• No grupo que votou em 2004 (\(primary2004=1\)) esse efeito será dado por \((\hat{\alpha} + \hat{\beta}_1 + \hat{\beta}_2 + \hat{\beta}_3) - (\hat{\alpha} + \hat{\beta}_1) = \hat{\beta}_2 + \hat{\beta}_3\). O primeiro termo ocorre quando as varfiáveis \(primary2004\) e \(messages\) são iguais a um e o segundo termo quando\(primary2004\) é igual a um e \(messages\) igual a zero.

• Para o grupo que não votou em 2004 (\(primary2004=0\)) o efeito médio do tratamento é dado por: \((\hat{\alpha} + \hat{\beta}_2) - \hat{\alpha} = \hat{\beta}_2\).
• A diferença do tratamento entre os que participaram e não participaram das primárias de 2004 é dada por: \(\hat{\beta}_2 + \hat{\beta}_3 - \hat{\beta}_2 = \hat{\beta}_3\)

• Até agora trabalhamos com variáveis categóricas, mas podemos analisar regressões com variáveis explicativas contínuas.
• Nesses casos é bom lembrar que a hipótese de linearidade pede que o impacto em \(Y\) do aumento de \(X\) em uma unidade não depende do valor de \(X\).
• Consideremos a idade como variável explicativa.

• Criar uma variável para idade.

social.neighbor <- social.neighbor %>%
  mutate(age = 2008 - yearofbirth)

summary(social.neighbor$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.00   43.00   52.00   51.82   61.00  108.00

• Para examinar como o efeito do tartamento, receber carta do tipo Neighbor muda com a idade vamos estimar uma regressão do tipo: \[Y=\alpha + \beta_1 \mbox{age} + \beta_2 \mbox{Neighbors} + \beta_3 \mbox{age}\times\mbox{Neighbors} + \varepsilon\]

• Usemos a mesma estratégia do exemplo anterior. O efeito da mensagem Neighbor em eleitores com \(x\) anos será: \((\hat{\alpha} + \hat{\beta}_1 x + \hat{\beta}_2 + \hat{\beta}_3 x) - (\hat{\alpha} + \hat{\beta}_1 x) = \hat{\beta}_2 + \hat{\beta}_3 x\).
• Em eleitores com idade \(x+1\) o efeito será: \((\hat{\alpha} + \hat{\beta}_1 (x+1) + \hat{\beta}_2 + \hat{\beta}_3 (x+1)) - (\hat{\alpha} + \hat{\beta}_1 (x+1)) = \hat{\beta}_2 + \hat{\beta}_3 (x+1)\).
• A diferença do tartamento em eleitores com \(x\) e \(x+1\) anos será: \((\hat{\beta}_2 + \hat{\beta}_3 (x+1)) - (\hat{\beta}_2 + \hat{\beta}_3 x) = \hat{\beta}_3\).

fit.age <- lm(primary2006 ~ age*messages, data = social.neighbor)
fit.age

## 
## Call:
## lm(formula = primary2006 ~ age * messages, data = social.neighbor)
## 
## Coefficients:
##           (Intercept)                    age      messagesNeighbors  
##             0.0894768              0.0039982              0.0485728  
## age:messagesNeighbors  
##             0.0006283

• O resultado sugere que a diferença no efeito do tratamento entre grupos de eleitores com diferença de um ano de idade é de 0,06%.
• Podemos usar o modelo para estimar o efeito para diferentes idades, por exemplo, 20, 40, 60 e 80 anos.
• Antes vamos ilustrar a função crossing() do pacote tidyr que será utilizada para criar as combinações entre idade e tartamento recebido.

crossing(age = seq(from=20, to=80, by=20), messages = c("Neighbors", "Control"))

## # A tibble: 8 × 2
##     age messages 
##   <dbl> <chr>    
## 1    20 Control  
## 2    20 Neighbors
## 3    40 Control  
## 4    40 Neighbors
## 5    60 Control  
## 6    60 Neighbors
## 7    80 Control  
## 8    80 Neighbors

ate.age <- tidyr::crossing(age = seq(from=20, to=80, by=20), 
                    messages = c("Neighbors", "Control")) %>%
  add_predictions(fit.age) %>%
  pivot_wider(names_from = messages, values_from = pred) %>%
  mutate(diff = Neighbors - Control)

ate.age

## # A tibble: 4 × 4
##     age Control Neighbors   diff
##   <dbl>   <dbl>     <dbl>  <dbl>
## 1    20   0.169     0.231 0.0611
## 2    40   0.249     0.323 0.0737
## 3    60   0.329     0.416 0.0863
## 4    80   0.409     0.508 0.0988

• Pesquisadores encontraram que a hipótese de linearidade não é adequada para modelar comparecimento eleitoral. A medida que a pessoa envelhece aumenta disposição a votar, mas, a partir do 60 ou 70 anos, cai a probabilidade de ir votar.
• Uma estratégia popular para enfrentar questões do tipo é modelar a variável de interesse, no caso comparecimento eleitoral, como função quadrática da variável explicativa, que no exemplo é a idade.

• Considere o modelo: \[Y=\alpha + \beta_2 \mbox{age}^2 + \beta_3 \mbox{Neighbors} + \beta_4 \mbox{age} \times \mbox{Neighbors} + \\ + \beta_5 \mbox{age}^2 \times \mbox{Neighbors} + \varepsilon\]

• Para colocar funções como quadrado na fórmula de uma regressão no R use a função I().

fit.age2 <- lm(primary2006 ~ age + I(age^2) + messages + age:messages + I(age^2):messages,
               data = social.neighbor)

fit.age2

## 
## Call:
## lm(formula = primary2006 ~ age + I(age^2) + messages + age:messages + 
##     I(age^2):messages, data = social.neighbor)
## 
## Coefficients:
##                (Intercept)                         age  
##                 -9.700e-02                   1.172e-02  
##                   I(age^2)           messagesNeighbors  
##                 -7.389e-05                  -5.275e-02  
##      age:messagesNeighbors  I(age^2):messagesNeighbors  
##                  4.804e-03                  -3.961e-05

• Em modelos complexos como o em fit.age2 interpretar os coeficientes é uma tarefa mais complicada, acaba sendo melhor usar o modelo para prever os valor da variável de resultado em diferentes cenários.
• Vamos usar a função data_grid() do pacote modelr para criar todas as combinações possíveis de idade e tratamento, estas combinações serão os cenários para os quais faremos as previsões.

social.neighbor %>%
  data_grid(age, messages)

## # A tibble: 172 × 2
##      age messages 
##    <dbl> <chr>    
##  1    22 Control  
##  2    22 Neighbors
##  3    23 Control  
##  4    23 Neighbors
##  5    24 Control  
##  6    24 Neighbors
##  7    25 Control  
##  8    25 Neighbors
##  9    26 Control  
## 10    26 Neighbors
## # ℹ 162 more rows

y.hat <- social.neighbor %>%
  data_grid(age, messages) %>%
  add_predictions(fit.age2)

head(y.hat)

## # A tibble: 6 × 3
##     age messages   pred
##   <dbl> <chr>     <dbl>
## 1    22 Control   0.125
## 2    22 Neighbors 0.159
## 3    23 Control   0.134
## 4    23 Neighbors 0.170
## 5    24 Control   0.142
## 6    24 Neighbors 0.182

• Para facilitar a apresentação dos resultados vamos fazer um gráfico com os valores previstos.

ggplot(y.hat, aes(age, pred)) +
  geom_line(aes(linetype = messages, color = messages)) +
  labs(color = "", linetype="",
       y = "Predicited turnout rate", x="Age") +
  xlim(20,90) +
  theme_classic()

• Também podemos fazer previsões para o efeito médio do tratamento por idade.

social.neighbor.ate <- y.hat %>%
  pivot_wider(names_from = messages, values_from = pred) %>%
  mutate(ate = Neighbors - Control)

head(social.neighbor.ate)

## # A tibble: 6 × 4
##     age Control Neighbors    ate
##   <dbl>   <dbl>     <dbl>  <dbl>
## 1    22   0.125     0.159 0.0338
## 2    23   0.134     0.170 0.0368
## 3    24   0.142     0.182 0.0397
## 4    25   0.150     0.192 0.0426
## 5    26   0.158     0.203 0.0454
## 6    27   0.166     0.214 0.0481

• Também podemos fazer previsões para o efeito médio do tratamento por idade.

head(social.neighbor.ate)

## # A tibble: 6 × 4
##     age Control Neighbors    ate
##   <dbl>   <dbl>     <dbl>  <dbl>
## 1    22   0.125     0.159 0.0338
## 2    23   0.134     0.170 0.0368
## 3    24   0.142     0.182 0.0397
## 4    25   0.150     0.192 0.0426
## 5    26   0.158     0.203 0.0454
## 6    27   0.166     0.214 0.0481

social.neighbor.ate %>%
  ggplot(aes(age, ate)) +
  geom_line() +
  labs(y="Estimated average treatment effect",
       x = "Age") +
  xlim(20,90) +
  ylim(0,0.1) +
  theme_classic()

• Na unidade 2 vimos que a associação entre variável de resultado e tratamento pode ser interpretada como causal na ausência de variáveis confundidoras. Esse costuma ser o caso em dados obtidos por experimentos aleatórios.
• Em estudos com dados observados a distribuição do tratamento não é aleatória, como resultado, confundidores, em vez do tratamento, podem explicar as diferenças observadas nos resultados dos grupos de tratamento e controle.
• Esse problema motiva a criação de uma série de estratégias para obter relações causais a partir de dados observados.

• Uma dessas técnica é a regressão com descontinuidade (RDD, do inglês: regression descontinuity design).
• Encerramos esta unidade apresentando um exemplo de uso de RDD.

• O objetivo é saber quanto um político pode aumentar a riqueza pessoal por estar ocupando um cargo.
• Os dados são referentes a membros do parlamento britânico. Os autores do estudo original coletaram informações sobre riqueza pessoal no momento da morte referentes a centenas de candidatos que disputaram eleições para o parlamento britânico entre 1950 e 1970. Os dados estão no data.frame MPs.
• O exemplo tem por base o artigo: Andrew C. Eggers e Jens Mainmueller (2009). MPs for sale? Returns to office in postwar British politics. American Political Science Review, v.103, n.4.

MPs <- read.csv("MPs.csv")

• Variáveis:
  • surname: sobrenome do candidato
  • firstname: primeiro nome do candidato
  • party: partido do candidato (labour ou tory)
  • ln.gross: log da riqueza bruta no momento da morte
  • ln.net: log da riqueza líquida no momento da morte
  • yob: ano de nascimento do candidato
  • yod: ano da morte do candidato
  • margin.pre: margem do candidato do partido na eleição anterior
  • region: região eleitoral
  • margin: margem de vitória (percentual dos votos)

• Simplesmente comparar a riqueza dos que foram eleitos com os que não foram eleitos não seria adequado porque os eleitos podem diferir dos não eleitos em várias características observáveis e não observáveis.
• No lugar dessa comparação, o RDD compara os candidatos que venceram por pouco com os que perderam por pouco.

• A ideia é avaliar se uma pequena mudança na margem de votos que leve de positiva para negativa (vitória para derrota) está associada a um grande salto positivo na riqueza do vencedor.
• Supondo que nada mais ocorre no ponto de descontinuidade, podemos identificar o efeito causal na riqueza de ser um parlamentar comparando candidatos que perderam por pouco com os que ganharam por pouco.
• A regressão é usada para prever a riqueza pessoal no ponto de descontinuidade.

• Um diagrama de dispersão com as retas de regressão é a melhor maneira de entender RDD.
• Vamos fazer o diagrama da variável de resultado (log da riqueza líquida ao morrer) contra a margem de vitória. No gráfico será colocada uma reta de regressão estimada com dados de candidatos vitoriosos (margem positiva) e outra estimada com dados de candidatos derrotados (margem negativa).
• Será feito um gráfico para cada partido.

• Começamos criando quatro subdivisões dos dados com as combinações partido (tory e labour) e margem (positiva e negativa). Depois faremos as regressões em cada conjunto.

labour_winners <- filter(MPs, party=="labour", margin > 0)
labour_losers <- filter(MPs, party=="labour", margin < 0)
tory_winners <- filter(MPs, party=="tory", margin > 0)
tory_losers <- filter(MPs, party=="tory", margin < 0)

labour_fit_win <- lm(ln.net ~ margin, data = labour_winners)
labour_fit_lose <- lm(ln.net ~ margin, data = labour_losers)
tory_fit_win <- lm(ln.net ~ margin, data = tory_winners)
tory_fit_lose <- lm(ln.net ~ margin, data = tory_losers)

• Na sequência adicionamos os valores previstos à cada conjunto.

y1_labour_win <- labour_winners %>%
  data_grid(margin) %>%
  add_predictions(labour_fit_win)

y2_labour_lose <- labour_losers %>%
  data_grid(margin) %>%
  add_predictions(labour_fit_lose)

y1_tory_win <- tory_winners %>%
  data_grid(margin) %>%
  add_predictions(tory_fit_win)

y2_tory_lose <- tory_losers %>%
  data_grid(margin) %>%
  add_predictions(tory_fit_lose)

head(y1_labour_win)

## # A tibble: 6 × 2
##    margin  pred
##     <dbl> <dbl>
## 1 0.00243  11.9
## 2 0.00372  11.9
## 3 0.00647  11.9
## 4 0.00938  12.0
## 5 0.0150   12.0
## 6 0.0161   12.0

ggplot() +
  geom_point(data = labour_winners, aes(margin, ln.net), shape = 1) +
  geom_point(data = labour_losers, aes(margin, ln.net), shape = 1) +
  geom_line(data = y1_labour_win, aes(margin, pred), color="blue", linewidth=1) +
  geom_line(data = y2_labour_lose, aes(margin, pred), color="blue", linewidth=1) +
  geom_vline(xintercept=0, linetype="dashed") +
  labs(title = "Labour",
       x = "Margin of victory",
       y = "log net wealth at death") +
  xlim(-0.5, 0.5) +
  ylim(6, 18) +
  theme_classic()

ggplot() +
  geom_point(data = tory_winners, aes(margin, ln.net), shape = 1) +
  geom_point(data = tory_losers, aes(margin, ln.net), shape = 1) +
  geom_line(data = y1_tory_win, aes(margin, pred), color="blue", linewidth=1) +
  geom_line(data = y2_tory_lose, aes(margin, pred), color="blue", linewidth=1) +
  geom_vline(xintercept=0, linetype="dashed") +
  labs(title = "Tory",
       x = "Margin of victory",
       y = "log net wealth at death") +
  xlim(-0.5, 0.5) +
  ylim(6, 18) +
  theme_classic()

• O resultado sugere que ser membro do parlamento aumenta a riqueza dos políticos Tory, mas não dos Labour.
• Podemos usar o R para calcular o tamanho do efeito. Para fazer isso vamos prever a riqueza quando a margem é zero usando os modelos para perdedor e vencedor, em seguida aplicamos exponencial, exp(), para reverter o logaritmo e finalmente calculamos a diferencça.
• Para adicionar a previsão dos dois modelo usamos a função spread_predictions() do pacote modelr.

spread_predictions(tibble(margin=0),
                   tory_fit_win, tory_fit_lose)

## # A tibble: 1 × 3
##   margin tory_fit_win tory_fit_lose
##    <dbl>        <dbl>         <dbl>
## 1      0         13.2          12.5

spread_predictions(tibble(margin=0),
                   tory_fit_win, tory_fit_lose) %>%
  mutate(rd_est = exp(tory_fit_win) - exp(tory_fit_lose)) %>%
  select(rd_est) %>%
  pull()

##        1 
## 255050.9

• O efeito estimado pelo modelo é que, para um Tory, ser membro do parlamento aumenta a riqueza em cerca de 250,000 pounds. Como a riqueza média de um Tory que não foi membro do parlamento é de 270,000 pounds o efeito é signficativo, ser membro do parlamento quase dobra a riqueza de um Tory.

• Para avaliar a validade interna da RDD podemos fazer testes de placebo.
• O teste consiste em encontrar situações onde teoricamente o efeito é zero e estimar o efeito para checar se é mesmo próximo de zero.
• No exemplo vamos estimar o efeito na margem de vitória do mesmo partido na eleição anterior. Como ser eleito no futuro não afeta o resultado da eleição passada o efeito deve ser zero.

• Primeiro fazemos as regressões:

tory_fit_win_placebo <- lm(margin.pre ~ margin, data = tory_winners)
tory_fit_lose_placebo <- lm(margin.pre ~ margin, data = tory_losers)

• O efeito estimado do tratamento é a diferença entre os dois interceptos. Para obter os interceptos vamos usar a função tidy() do pacote tidyverse.

win_intercept <- tidy(tory_fit_win_placebo) %>%
  filter(term == "(Intercept)") %>%
  select(estimate) %>%
  pull()

lose_intercept <- tidy(tory_fit_lose_placebo) %>%
  filter(term == "(Intercept)") %>%
  select(estimate) %>%
  pull()

win_intercept - lose_intercept

## [1] -0.01725578

• Como esperado, o efeito foi bem próximo de zero.

• RDD pode ajudar a resolver problemas enfrentados por pesquisas usando dados observados, porém, a melhora na validade interna pode vir ao custo de piora na validade externa. Isso acontece porque os efeitos causais estimados só se aplicam na fronteira da discontinuidade.
• No nosso exemplo, é possível que o efeito damargem sobre a riqueza seja diferente para candidatos que ganham as eleições por uma larga margem.

• Regressão com descontinuidade é uma estratégia de pesquisa para realizar inferência causal em estudos com dados observados na presença de possíveis confundidores. Um RDD supõe que a mudança do resultado no ponto de descontinuidade pode ser atribuída a variação na variável de tratamento. Resultados obtidos com RDD ganham em validade interna, mas perdem em validade externa.