• Na sequência dos atentados de 11 de Setembro de 2001 os Estados Unidos invadiram o Afeganistão com o objetivo de destruir a al-Qaeda. Em 2003 a OTAN entrou na guerra e mandou tropas com o nome de Força Internacional de Apoio à Segurança (ISAF, do inglês International Security Assistance Force).
• A ISAF adotou uma campanha para conquistar corações e mentes da população, a campanha combinava assistência econômica, prestação de serviços e proteção aos civis.
• Para avaliar a campanha é essencial entender e medir as experiências e sentimentos dos civis durante o período de guerra. Não bastasse as dificuldades inatas de medir sentimentos, as condições de guerra podiam tornar perigoso revelar os sentimentos e as experiencias respondendo as questões das pesquisas.

• Um grupo de cientistas sociais resolveu fazer uma pesquisa de opinião no sul do Afeganistão, o coração da guerra de insurgência.
• São dois os artigos que fundamentam esse exemplo:
  • Jason Lyall, Graemer Blair e Kosuke Imai. Explaining support for combatants during wartime: a survey experiment in Afghanistan. American Political Science Review, v.107, 2013.
  • Graeme Blair, Kosuke Imai e Jason Lyall. Comparing and combining list and endorsement experiments: evidence from Afghanistan. American Journal of Political Science, v.58, 2014.

• A pesquisa foi aplicada em uma amostra de 2754 cidadãos do Afeganistão entre janeiro e fevereiro de 2011.
• A taxa de participação foi de 89%, ou seja, foram contactados 3097 homens, mas 343 se recusaram a responder.
• Por questões locais, apenas homens responderam às pesquisas.

• Variáveis:
  • province: província onde vive quem respondeu.
  • district: distrito onde vive quem respondeu.
  • village.id: o ID da vila onde vive quem respondeu.
  • age: idade.
  • educ.years: anos de educação.

• Variáveis:
  • employed: se quem respondeu está empregado.
  • income: renda mensal (cinco níveis)
  • violent.ex.ISAF: se experimentou violência pelo ISAF
  • violent.ex.taliban: se experimentou violência pelo taliban
  • list.group: grupo designado aleatoriamente para list experiment (control, ISAF, taliban)
  • list.response: repostas para as questões do list experiment (0-4)

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

## 'data.frame':    2754 obs. of  11 variables:
##  $ province           : chr  "Logar" "Logar" "Logar" "Logar" ...
##  $ district           : chr  "Baraki Barak" "Baraki Barak" "Baraki Barak" "Baraki Barak" ...
##  $ village.id         : int  80 80 80 80 80 80 80 80 80 35 ...
##  $ age                : int  26 49 60 34 21 18 42 39 20 18 ...
##  $ educ.years         : int  10 3 0 14 12 10 6 12 5 10 ...
##  $ employed           : int  0 1 1 1 1 1 1 1 1 0 ...
##  $ income             : chr  "2,001-10,000" "2,001-10,000" "2,001-10,000" "2,001-10,000" ...
##  $ violent.exp.ISAF   : int  0 0 1 0 0 0 0 0 0 0 ...
##  $ violent.exp.taliban: int  0 0 0 0 0 0 0 1 0 0 ...
##  $ list.group         : chr  "control" "control" "control" "ISAF" ...
##  $ list.response      : int  0 1 1 3 3 2 1 3 1 2 ...

head(afghan)

##   province     district village.id age educ.years employed       income
## 1    Logar Baraki Barak         80  26         10        0 2,001-10,000
## 2    Logar Baraki Barak         80  49          3        1 2,001-10,000
## 3    Logar Baraki Barak         80  60          0        1 2,001-10,000
## 4    Logar Baraki Barak         80  34         14        1 2,001-10,000
## 5    Logar Baraki Barak         80  21         12        1 2,001-10,000
## 6    Logar Baraki Barak         80  18         10        1         <NA>
##   violent.exp.ISAF violent.exp.taliban list.group list.response
## 1                0                   0    control             0
## 2                0                   0    control             1
## 3                1                   0    control             1
## 4                0                   0       ISAF             3
## 5                0                   0       ISAF             3
## 6                0                   0       ISAF             2

library(tidyverse)
afghan %>%
  select(age, educ.years, employed, income) %>%
  summary()

##       age          educ.years        employed         income         
##  Min.   :15.00   Min.   : 0.000   Min.   :0.0000   Length:2754       
##  1st Qu.:22.00   1st Qu.: 0.000   1st Qu.:0.0000   Class :character  
##  Median :30.00   Median : 1.000   Median :1.0000   Mode  :character  
##  Mean   :32.39   Mean   : 4.002   Mean   :0.5828                     
##  3rd Qu.:40.00   3rd Qu.: 8.000   3rd Qu.:1.0000                     
##  Max.   :80.00   Max.   :18.000   Max.   :1.0000

• Eu definiria income como fator…

afghan %>%
  select(age, educ.years, employed, income) %>%
  mutate(income = as.factor(income)) %>%
  summary()

##       age          educ.years        employed                  income    
##  Min.   :15.00   Min.   : 0.000   Min.   :0.0000   10,001-20,000  : 616  
##  1st Qu.:22.00   1st Qu.: 0.000   1st Qu.:0.0000   2,001-10,000   :1420  
##  Median :30.00   Median : 1.000   Median :1.0000   20,001-30,000  :  93  
##  Mean   :32.39   Mean   : 4.002   Mean   :0.5828   less than 2,000: 457  
##  3rd Qu.:40.00   3rd Qu.: 8.000   3rd Qu.:1.0000   over 30,000    :  14  
##  Max.   :80.00   Max.   :18.000   Max.   :1.0000   NA's           : 154

• … mas vou seguir o livro.

count(afghan, income)

##            income    n
## 1   10,001-20,000  616
## 2    2,001-10,000 1420
## 3   20,001-30,000   93
## 4 less than 2,000  457
## 5     over 30,000   14
## 6            <NA>  154

unique(afghan$income)

## [1] "2,001-10,000"    NA                "10,001-20,000"   "less than 2,000"
## [5] "20,001-30,000"   "over 30,000"

• A média de idade dos que responderam as questões é de 32 anos.
• A grande maioria tem pouca educação (a mediana é 1)
• Cerca de 60% estavam desempregados.
• A maioria dos que responderam tem renda mensal inferior a 10.000 afghani, que na época era equivalente a 200 dólares.

• Para avaliar como a violência atingiu a população civil, foi feita a seguinte pergunta: No último ano, você ou alguém de sua família foi atingido por conta da ação de forças entrangeiras ou do taliban?
• Pergunta orginal: Over the past year, have you or anyone in your family suffered harm due to the actions of the Foreign Forces/the Taliban?

• Na pesquisa foi registrado que harm, que traduzi como atingido, se refere a injúrias físicas bem como danos ao patrimônio.
• Vamos calcular a proporção dos que foram atingidos por ações da ISAF e/ou do Taliban.

harm_props <- afghan %>%
  group_by(violent.exp.ISAF, violent.exp.taliban) %>%
  count() %>%
  ungroup() %>%
  mutate(prop = n/sum(n))

## # A tibble: 9 × 4
##   violent.exp.ISAF violent.exp.taliban     n    prop
##              <int>               <int> <int>   <dbl>
## 1                0                   0  1330 0.483  
## 2                0                   1   354 0.129  
## 3                0                  NA    22 0.00799
## 4                1                   0   475 0.172  
## 5                1                   1   526 0.191  
## 6                1                  NA    22 0.00799
## 7               NA                   0     7 0.00254
## 8               NA                   1     8 0.00290
## 9               NA                  NA    10 0.00363

• Sem usar ungroup() a soma da última linha seria feita dentro de cada grupo.

afghan %>%
  group_by(violent.exp.ISAF, violent.exp.taliban) %>%
  count() %>%
#  ungroup() %>%
  mutate(prop = n/sum(n))

## # A tibble: 9 × 4
## # Groups:   violent.exp.ISAF, violent.exp.taliban [9]
##   violent.exp.ISAF violent.exp.taliban     n  prop
##              <int>               <int> <int> <dbl>
## 1                0                   0  1330     1
## 2                0                   1   354     1
## 3                0                  NA    22     1
## 4                1                   0   475     1
## 5                1                   1   526     1
## 6                1                  NA    22     1
## 7               NA                   0     7     1
## 8               NA                   1     8     1
## 9               NA                  NA    10     1

• Proporção atingida pela ISAF.

ISAF_harm_prop <- harm_props %>%
  filter(violent.exp.ISAF == 1) %>%
  summarize(harm_prop = sum(prop)) %>%
  pull()

ISAF_harm_prop

## [1] 0.3714597

• Proporção atingida pelo taliban

talib_harm_prop <- harm_props %>%
  filter(violent.exp.taliban == 1) %>%
  summarize(harm_prop = sum(prop)) %>%
  pull()

talib_harm_prop

## [1] 0.3224401

• Proporção atingida pelo taliban e pela ISAF

both_harm_prop <- harm_props %>%
  filter(violent.exp.taliban == 1 & violent.exp.ISAF == 1) %>%
  summarize(harm_prop = sum(prop)) %>%
  pull()

both_harm_prop

## [1] 0.1909949

• Em muitos questionários aparecem perguntas não respondidas, mesmo em bases com dados observados é comum que alguns dados estejam faltando.
• Quando o dado está faltando o R registra como NA.

afghan %>%
  select(income) %>%
  slice(1:10)

##           income
## 1   2,001-10,000
## 2   2,001-10,000
## 3   2,001-10,000
## 4   2,001-10,000
## 5   2,001-10,000
## 6           <NA>
## 7  10,001-20,000
## 8   2,001-10,000
## 9   2,001-10,000
## 10          <NA>

• É fundamental ter certeza que todos os dados ausentes são reconhecidos pelo R como ausentes e registrados como NA. Algumas vezes garantir que isso aconteça é trabalhoso, algumas bases registram NAs como espaços em branco, -, –, ou números estranhos como -999. Um NA não reconhecido pelo R pode comprometer toda a análise de dados.
• A função is.na() identifica os NAs presentes nos dados.

afghan %>%
  select(income) %>%
  slice(1:10) %>%
  is.na()

##       income
##  [1,]  FALSE
##  [2,]  FALSE
##  [3,]  FALSE
##  [4,]  FALSE
##  [5,]  FALSE
##  [6,]   TRUE
##  [7,]  FALSE
##  [8,]  FALSE
##  [9,]  FALSE
## [10,]   TRUE

• Como FALSE equivale a 0 e TRUE equivale a 1, a soma dos resultados da função is.na() é o número de NAs e a média é a proporção de NAs.

afghan %>%
  summarize(n_missing = sum(is.na(income)),
            p_missing = mean(is.na(income)))

##   n_missing  p_missing
## 1       154 0.05591866

• No R os NAs são contagiosos, ou seja, operações com NAs retornam NAs.

5 + NA

## [1] NA

5*NA

## [1] NA

• Por serem contagiosos o pesquisador é forçado a decidir o que fazer com os NAs. Em alguns casos ignorá-los é seguro, em outros não. Cabe ao pesquisador decidir o que fazer.
• Por exemplo, suponha dados de estações que medem a velocidade do vento e que os aparelhos não conseguem medir velocidades muito altas. Ignorar os NAs vai levar a subestimar a velocidade do vento na região.

• Muita funções do R possuem o argumento na.rm que permite ignorar NAs. Essa não é uma prática errada, mas o pesquisador deve estar ciente dos riscos ao escolher ignorar os NAs.

x <- c(1,2,3,4,NA)
sum(x)

## [1] NA

sum(x, na.rm = TRUE)

## [1] 10

x <- c(1,2,3,4,NA)
mean(x)

## [1] NA

mean(x, na.rm = TRUE)

## [1] 2.5

• A função filter() combinada com a função is.na() oeferece uma maneira rápida de tirar NAs de uma análise.

afghan %>%
  filter(!is.na(violent.exp.ISAF), !is.na(violent.exp.taliban)) %>%
  group_by(violent.exp.ISAF, violent.exp.taliban) %>%
  count() %>%
  ungroup() %>%
  mutate(prop = n/sum(n)) %>%
  arrange(prop)

## # A tibble: 4 × 4
##   violent.exp.ISAF violent.exp.taliban     n  prop
##              <int>               <int> <int> <dbl>
## 1                0                   1   354 0.132
## 2                1                   0   475 0.177
## 3                1                   1   526 0.196
## 4                0                   0  1330 0.495

• A função arrange() organizou os dados da menor para maior proporção, podemos organizar em ordem decrescente usando a função desc().

afghan %>%
  filter(!is.na(violent.exp.ISAF), !is.na(violent.exp.taliban)) %>%
  group_by(violent.exp.ISAF, violent.exp.taliban) %>%
  count() %>%
  ungroup() %>%
  mutate(prop = n/sum(n)) %>%
  arrange(desc(prop))

## # A tibble: 4 × 4
##   violent.exp.ISAF violent.exp.taliban     n  prop
##              <int>               <int> <int> <dbl>
## 1                0                   0  1330 0.495
## 2                1                   1   526 0.196
## 3                1                   0   475 0.177
## 4                0                   1   354 0.132

• As funções filter() e is.na() também poder ser usadas para identificar NAs em mais de uma variável, por exemplo, podemos encontrar quantos NAs para quem sofreu violência com ISAF ou taliban.

missing_prop <- harm_props %>%
  filter(is.na(violent.exp.ISAF) | is.na(violent.exp.taliban)) %>%
  summarize(missing_prop = sum(prop)) %>%
  pull

missing_prop

## [1] 0.02505447

• As funções na.omit(), do R básico, e drop_na(), do tidyverse excluem as observações onde aparece NA em pelo menos uma variável. Repare que as funções excluem toda a linha mesmo que esteja faltando dados para apenas uma variável.

nrow(afghan)

## [1] 2754

afghan.sub <- na.omit(afghan)
nrow(afghan.sub)

## [1] 2554

nrow(afghan)

## [1] 2754

afghan.sub2 <- drop_na(afghan)
nrow(afghan.sub2)

## [1] 2554

afghan %>% drop_na(income) %>% nrow()

## [1] 2600

• Além de estatísitcas descritivas podemos passar informações sobre a distribuição de uma variável de interesse usando gráficos.
• Na introdução vimos como usar o ggplot2 para criar gráficos por etapas, agora vamos usar esse conhecimento para passar informações a respeito de como os civis afegãos foram atingidos pela ISAF e pelo taliban.

• Para começar vamos ilustrar as respostas sobre para a pergunta sobre ter sido atingido pela ISAF, queremos comparar as repostas afirmativas, negativas e os que não responderam.
• Gráficos de barras fazem bem esse trabalho (muito cuidado com gráficos de pizza!). Na introdução usamos a função geom_col() para fazer gráficos de barras, aqui vamos usar a função geom_bar(). A função geom_col() é uma versão simplificada da geom_bar().

• Vítimas da ISAF

ggplot(data = afghan, aes(x=as.factor(violent.exp.ISAF))) +
  geom_bar(aes(y= after_stat(prop))) +
  scale_x_discrete(labels = c("No harm", "Harm", "Nonresponse")) +
  ylab("Proportion of respondents") +
  xlab("Response category") +
  ggtitle("Civilian Victimization by the ISAF")

• No livro a função aes() no arumento da função geom_bar() é chamada da seguinte forma: geom_bar(aes(y=..prop.., group=1)) o argumento ..prop.. foi substituido por after_stat(prop) para ficar compatível com as versões mais modernas do ggplot2. O argumento group = 1 é um “truque” para fazer o agrupamento pelos níveis de violent.exp.ISAF. Por fim, note que os autores transformaram a variável violent.exp.ISAF em fator.
• As barras podem ser feitas por vários critérios, o argumento prop na função after_stat() diz para usar proporções. Para ilustrar a função after_stat() vou refazer o gráfico usando o argumento count.

ggplot(data = afghan, aes(x=as.factor(violent.exp.ISAF))) +
  geom_bar(aes(y= after_stat(count))) +
  scale_x_discrete(labels = c("No harm", "Harm", "Nonresponse")) +
  ylab("Proportion of respondents") +
  xlab("Response category") +
  ggtitle("Civilian Victimization by the ISAF")

• Vítimas do taliban

ggplot(data = afghan, aes(x=as.factor(violent.exp.taliban))) +
  geom_bar(aes(y= after_stat(prop))) +
  scale_x_discrete(labels = c("No harm", "Harm", "Nonresponse")) +
  ylab("Proportion of respondents") +
  xlab("Response category") +
  ggtitle("Civilian Victimization by the ISAF")

• Os gráficos ilustram bem cada uma das distribuições, por exemplo, é fácil ver que a proporção de No harm é maior do que de haarm nas duas distribuições.
• Porém, não é fácil comparar um distribuição com a outra por meio dos gráficos. Por exemplo, fica difícil saber se houve mais harm para ISAF ou taliban apenas olhando os gráficos.
• Uma nova arrumação dos dados permite criar um gráfico onde a comparação é mais fácil.

afghan %>%
  pivot_longer(violent.exp.ISAF:violent.exp.taliban, 
               names_to ="harming group", values_to="harm") %>%
  ggplot(aes(x=as.factor(harm))) +
  geom_bar(aes(y=after_stat(prop), fill=harming_group, group=harming_group),
           position = "dodge") +
  scale_x_discrete(labels=c("No harm", "Harm", "Nonresponse")) +
  scale_fill_discrete(name="Harming group", labels=c("ISAF", "Taliban")) +
  labs(title = "Civilian victimization",
       x = "Response category",
       y = "Proportion of respondents")

• Para entender melhor o que aconteceu é válido observar o data.frame original e como foi modificado nas linhas que antecedem a função ggplot()

afghan %>%
  pivot_longer(violent.exp.ISAF:violent.exp.taliban, 
               names_to ="harming group", values_to="harm")

• Que tal mudar a paleta de cores e o tema?

afghan %>%
  pivot_longer(violent.exp.ISAF:violent.exp.taliban, 
               names_to ="harming_group", values_to="harm") %>%
  ggplot(aes(x=as.factor(harm))) +
  geom_bar(aes(y=after_stat(prop), fill=harming_group, group=harming_group),
           position = "dodge") +
  scale_x_discrete(labels=c("No harm", "Harm", "Nonresponse")) +
  scale_fill_brewer(name="Harming group", labels=c("ISAF", "Taliban"), 
                      palette = "Blues") +
  labs(title = "Civilian victimization",
       x = "Response category",
       y = "Proportion of respondents") +
  theme_classic()

• O histograma é um método popular para ilustrar distribuições de variáveis numéricas. Suponha, por exemplo, que queiramos dar uma olhada na idade de quem respondeu às enquetes. Como temos muitos valores um gráfico de barras ficaria confuso e não nos daria muita informação.

• Uma alternativa é dividir as idades em intervalos iguais (bins), então contamos quantas observações estão em cada intervalo (janela) e calculamos a densidade de acordo com a fórmula: \[\mbox{densidade}=\frac{\mbox{proporção de observações na janela}}{\mbox{tamanho da janela}}\]
• Em geral estamos menos preocupados com o valor exato da densidade do que em comparar a densidade das várias janelas.
• A função geom_histogram() do ggplot2 faz o histograma.

afghan %>%
  ggplot(aes(x = age)) +
  geom_histogram(aes(y= after_stat(density)), binwidth=5, boundary=0) +
  scale_x_continuous(breaks = seq(20,80,by=10)) +
  labs(title = "Distribution of respondent's age",
       y = "Density",
       x = "Age") +
  theme_classic()

• Em um histograma a área do retângulo associado a cada janela representa a proporção de observações na janela.
• Desta forma, em geral, interpretamos a escala da densidade, a unidade no eixo vertical, como percentagem pela unidade do eixo horizontal. No exemplo, a densidade é medida como percentual por idade.
• Por consequência, as alturas dos retângulos associados às janelas pode somar mais do que um, já as áreas dos retângulos somam um.

• Um histograma divide os dados em janelas onde a área de cada janela representa a proporção de observações que pertencem à janela. A altura de cada janela é chamada de densidade, que é igual a proporção de observações em cada janela dividida pela largura da janela. O histograma aproxima a distribuição de uma variável.

• Histograma como número de observações no lugar da densidade.

afghan %>%
  ggplot(aes(x = age)) +
  geom_histogram(aes(y= after_stat(count)), binwidth=5, boundary=0) +
  scale_x_continuous(breaks = seq(20,80,by=10)) +
  labs(title = "Distribuição das idades",
       y = "Número de observações",
       x = "Idade") +
  theme_classic()

• A função geom_vline() coloca linhas verticais em um gráfico (para linhas horizontais use a função geom_hline()) e a função annotate() faz anotações no gráfico.
• Vamos fazer o histograma dos anos de educação destacando a mediana.

afghan %>%
  ggplot(aes(x = educ.years, y= after_stat(density))) +
  geom_histogram(binwidth=1, center=0) +
  geom_vline(xintercept = median(afghan$educ.years)) +
  annotate(geom="text", x = median(afghan$educ.years),
           y=0.4, label="Median", hjust=-0.1)
  labs(title = "Distribution of respondent's education",
       y = "denssity",
       x = "Years of education") +
  theme_classic()

• Outro tipo de figura muito usado para ilustrar distribuições é o boxplot, diagrama de caixas. O boxplot é particularmente útil para comparar distribuições de diferentes variáveis.
• O boxplot desenha um retângulo (caixa) com 50% das observações. A altura da caixa é o intervalo interquartílico, a mediana é marcada por uma linha na caixa, abaixo da caixa tem um linha que mostra 1,5 IQR abaixo do primeiro quartil e acima da caixa uma linha que mostra 1,5 acima do terceiro quartil. Os pontos abaixo e acima da linha, pontos extremos ou outliers, são destacado na figura.
• Na ausência de pontos extremos baixos, a linha inferior acaba na menor observação. Da mesma forma, na ausência de pontos extremos altos, a linha superior acaba na maior observação.

• Boxplot da idade

afghan %>%
  ggplot(aes(y=age)) +
  geom_boxplot() +
  labs(y="Age", x=NULL, title = "Distribution of age")

• Como foi dito, o boxplot é particularmente útil para comparar distribuições. Por exemplo, suponha que você quer comparar a distribuição dos anos de educação nas diversas províncias onde a pesquisa foi aplicada.

afghan %>%
  ggplot(aes(y=educ.years, x=province)) +
  geom_boxplot() +
  labs(y="Years of education", x = "Province", title = "Education by province")

• É fácil ver que os participantes das províncias de Helmand e Uruzgam possuem menos anos de educação do que os das outras províncias, talvez não por coincidência sejam os que mais reportam ter sido alvo da ISAF e do taliban.

afghan %>%
  group_by(province) %>%
  summarize(violence.exp.taliban = mean(violent.exp.taliban, na.rm=TRUE),
            violence.exp.ISAF = mean(violent.exp.ISAF, na.rm=TRUE))

## # A tibble: 5 × 3
##   province violence.exp.taliban violence.exp.ISAF
##   <chr>                   <dbl>             <dbl>
## 1 Helmand                0.504              0.541
## 2 Khost                  0.233              0.242
## 3 Kunar                  0.303              0.399
## 4 Logar                  0.0802             0.144
## 5 Uruzgan                0.455              0.496

• Um diagrama de caixa ou boxplotmostra a distribuição de uma variável indicando a mediana, primeiro e terceiro quartis, e os pontos extremos. Permite uma comparação entre as distribuições de um conjunto de variáveis.

• A maneira mais fácil de salvar figuras é usando os menus do RSTudio, mas, para quem não gosta de seguir menus de forma alguma, existe a função ggsave(). Deixo por conta de vocês olhar o help desta função.

• Algumas vezes, por questões de estética ou para facilitar comprações, desejamos colocar figuras uma ao lado da outra. A função grid.arrange() do pacote gridExtra facilita esse trabalho.

library(gridExtra)

age_hist <- afghan %>%
  ggplot(aes(x=age)) +
  geom_histogram(aes(y= after_stat(density)), binwidth=5, boundary=0) +
  scale_x_continuous(breaks = seq(20,80,by=10)) +
  labs(title = "Distribution of respondent's age",
       y = "Density",
       x = "Age") +
  theme_classic()

educ_hist <- afghan %>%
  ggplot(aes(x = educ.years, y= after_stat(density))) +
  geom_histogram(binwidth=1, center=0) +
  geom_vline(xintercept = median(afghan$educ.years)) +
  annotate(geom="text", x = median(afghan$educ.years),
           y=0.4, label="Median", hjust=-0.1) +
  labs(title = "Distribution of respondent's education",
       y = "denssity",
       x = "Years of education") +
  theme_classic()

grid.arrange(age_hist, educ_hist, ncol=2)

• Pesquisa por amostragem é uma das principais ferramentas para cientistas sociais que fazem trabalhos aplicados. Normalmente são usadas para avaliar comportamento ou opinião da população quando as informações não estão disponíveis em outras fontes como, por exemplo, registros administrativos.
• Na pesquisa por amostragem o pesquisador seleciona um subconjunto da população, chamado de amostra, para entender características da população como um todo.
• Pesquisa por amostragem é diferente de censo, esse último tem como objetivo abranger toda a população.

• O uso de amostras permite obter informações sobre uma grande população entrevistando uma pequena fração dessa população.
• No exemplo do Afeganistão, uma amostra de 2754 pessoas foi usada para inferir experiências e comportamentos de uma população de mais de 15 milhões de pessoas.

• Nos EUA, amostras com pouco mais de mil pessoas costumam ser usadas para inferir a opinião de mais de 200 milhões de adultos.
• A “mágica” da amostra deve muito ao uso da escolha aleatória, aqui vamos tratar da amostragem aleatória simples.

• Antes do uso de amostras aleatórias, pesquisadores costumavam construir amostras por cotas. Grosso modo, o método consiste em determinar algumas caracteríticas da amostra para reproduzir a da população. Por exemplo, se 20% da população tem nível superior, então 20% da amostra deveria ser composta por pessoas com nível superior. Para definição das cotas é comum utilizar de multiplas variáveis.

• O problema das cotas é semelhante ao enfrentado por estudos com dados observados. Mesmo que a amostra reproduza caracterísitcas observadas da população que forma usadas para definir as cotas, é possível que a amostra tenha diferenças relevantes da população em termos de características não observadas.
• O uso de amostras aleatórias elimita esse viés de seleção da amostra garantindo que a amostra selecionada seja representativa da população alvo.

• Amostragem aleatória simples (SRS, do inglês: simple random sample) é o método mais simples de amostragem probabilística que evita problemas de viés de seleção da amostra por meio de escolha aleatória de unidades da população. No SRS, o número predeterminado de unidades é selecionado aleatoriamente e sem reposição de uma população-alvo, cada unidade tem a mesma probabilidade de ser selecionada. A amostra resultante é representativa da população em termos de características observadas e não observadas.

• O método de amostragem aleatória simples muitas vezes não é viável por conta de questões logísticas.
• Imagine que você vai aplicar esse método para fazer uma pesquisa aqui no DF. Onde conseguir uma lista com todos os moradores? Você pode ter pensado na lista telefônica, mas e quem não tem telefone? E quem tem vários telefones?
• Outras vezes pode ser muito difícil ou muito caro aplicar o questionário em alguns dos indivíduos escolhidos aleatóriamente.
• Por conta disso existe métodos mais sofisitcados de amostragem aleatórias.

• Para pesquisa do Afeganistão foi usado método de amostragem por cluster multiestágios. Esse método consiste em primeiro selecionar unidades grandes e depois escolher aleatoriamente unidades menores dentro das unidades maiores.
• No caso de Afeganistão, para cada um das cinco províncias de interesse foram selecionadas amostras de distritos e depois de vilas. Em cada vila os entrevistados foram escolhidos de maneira aproximadamente aleatória (a escolha levou em conta a localização na vila).

• O arquivo afghan-village.csv contém dados de altitude e população das vilas do Afeganistão.
• Vamos comparar essas propriedade nas vilas que entraram e não entraram na amostra, mas antes, vamos carregar e dar uma olhada nos dados.

afghan.village <- read.csv("afghan-village.csv")
head(afghan.village)

##   altitude population village.surveyed
## 1  1959.08        197                1
## 2  2425.88        744                0
## 3  2236.60        179                1
## 4  1691.76        225                0
## 5  1928.04        379                0
## 6  1194.56        617                0

• Histograma da altitude

afghan.village %>%
  ggplot(aes(x=altitude)) +
  geom_histogram(aes(y=after_stat(density))) +
  labs(x = "Altitude (metros)",
       y= "densidade")

• Histograma da população

afghan.village %>%
  ggplot(aes(x=population)) +
  geom_histogram(aes(y=after_stat(density))) +
  labs(x = "População",
       y= "densidade")

• Histograma do log da população

afghan.village %>%
  ggplot(aes(x=log(population))) +
  geom_histogram(aes(y=after_stat(density))) +
  labs(x = "População (log)",
       y= "densidade")

• Comparação entre as altitudes das vilas na amostra e fora da amostra

afghan.village %>%
  ggplot(aes(x = as.factor(x=village.surveyed),
             y = altitude)) +
  geom_boxplot() +
  scale_x_discrete(labels = c("Fora da amostra", "Na amostra")) +
  labs(y = "Altitude (metros)", x=NULL) +
  theme_classic()

• Comparação entre as populações das vilas na amostra e fora da amostra

afghan.village %>%
  ggplot(aes(x = as.factor(x=village.surveyed),
             y = log(population))) +
  geom_boxplot() +
  scale_x_discrete(labels = c("Fora da amostra", "Na amostra")) +
  labs(y = "População (log)", x=NULL) +
  theme_classic()

• Muitas vezes o cientista social é desafiado a desenvolver medidas para descrever ou entender comportamentos, atitudes e caracterísitcas não observadas das pessoas. Por exemplo, quantificar a ideologia de atores políticos como legisladores e juízes a partir do comportamento observado desses atores.
• Como medir o quão liberal ou conservador (lembre que o livro texto é “americano”) é um político ou, talvez mais próximo da realidade brasileira, quão de esquerda ou de direita é um senador?

• Uma estratégia adotada por cientistas políticos é observar os votos do congressista.
• O modelo de “espaço de votações”, spatial voting model, relacicona ideologia e votos de congressistas. O modelo apresenta duas dimensões: liberalismo/conservadorismo econômico e racial (pesquisadores identificaram essas como as princípais características ideológicas nos EUA do pós-guerra).
• Pesquisas mostram que o liberalismo/conservadorismo econômico tem mais poder explicativo para o voto do que o racial.

• De acordo com o modelo o congressista tem um ponto ideal no plano definido pelas dimensões relevantes, no caso liberalismo/conservadorismo racial e econômico.
• O congressista está mais propenso a votar contra uma proposta se o ponto ideal dele for mais próximo do status quo do que da proposta, vale o contrário.
• Propostas controversas revelam muito sobre a ideologia do congressista e propostas unânimes não nos dão informações sobre a ideologia do congressista.

• O exemplo de “medir ideologia” será usado para ilustrar relações bivariadas.
• A referência para o exemplo é: Nolan McCarty, Keiht T. Poole e Howard Rosenthal (2006). Polarized America: The Dance of Ideology and Unequal Richies. MIT Press.

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

##   congress district   state    party        name dwnom1 dwnom2
## 1       80        0     USA Democrat      TRUMAN -0.276  0.016
## 2       80        1 ALABAMA Democrat  BOYKIN  F. -0.026  0.796
## 3       80        2 ALABAMA Democrat   GRANT  G. -0.042  0.999
## 4       80        3 ALABAMA Democrat ANDREWS  G. -0.008  1.005
## 5       80        4 ALABAMA Democrat   HOBBS  S. -0.082  1.066
## 6       80        5 ALABAMA Democrat   RAINS  A. -0.170  0.870

• Variáveis:
  • name: nome do congressista.
  • state: estado do congressista.
  • party: partido do congressista.
  • congress: número da legislatura.
  • dwnom1: DW-NOMINATE score (dimensão 1).
  • dwnom2: DW-NOMINATE score (dimensão 2).
• Os dados vão da 80o legislatura (1947-48) ate a 112o legislatura (2011-2012)

• Vamos começar fazendo um gráfico de dispersão, scatter plot, para a 80o e 112o legislaturas. Republicanos e democratas serão caracterizados pela cor e pelo formato dos pontos.

plot_80 <- congress %>%
  filter(congress == 80) %>%
  ggplot(aes(x=dwnom1, y=dwnom2)) +
  geom_point(aes(shape=party, color = party), show.legend = FALSE) +
  scale_color_manual(values = c(Democrat = "blue",
                                Republican = "red",
                                Other = "green")) +
  scale_shape_manual(values = c(Democrat = "square",
                                Republican = "triangle",
                                Other = "circle")) +
  annotate("text",x=-1, y=1, label="Democrats", color = "blue") +
  annotate("text",x=0.5, y=-1, label="Republicans", color = "red") +
  scale_y_continuous("Racial liberalism/conservatism", limits = c(-1.5, 1.5)) +
  scale_x_continuous("Economic liberalism/conservatism", limits = c(-1.5, 1.5)) +
  ggtitle("80th Congress") + 
  coord_fixed() +
  theme_classic()

plot_112 <- congress %>%
  filter(congress == 112) %>%
  ggplot(aes(x=dwnom1, y=dwnom2)) +
  geom_point(aes(shape=party, color = party), show.legend = FALSE) +
  scale_color_manual(values = c(Democrat = "blue",
                                Republican = "red",
                                Other = "green")) +
  scale_shape_manual(values = c(Democrat = "square",
                                Republican = "triangle",
                                Other = "circle")) +
  annotate("text",x=-1, y=1, label="Democrats", color = "blue") +
  annotate("text",x=0.5, y=-1, label="Republicans", color = "red") +
  scale_y_continuous("Racial liberalism/conservatism", limits = c(-1.5, 1.5)) +
  scale_x_continuous("Economic liberalism/conservatism", limits = c(-1.5, 1.5)) +
  ggtitle("80th Congress") + 
  coord_fixed() +
  theme_classic()

grid.arrange(plot_80, plot_112, ncol=2)

• Note que as diferenças entre os partidos nas questões raciais são bem mais importantes no 80o Congresso do que no 112o Congresso. Vale o contrário nas questões econômicas.

• O diagrama de dispersão compara duas variáveis medidas para o mesmo conjunto de unidades em um gráfico com o valor de uma variável em um eixo e o de outra variável no outro eixo.

• O próximo exercício será comparar a mediana da primeira dimensão do DW-NOMINATE (liberalismo/conservadorismo econômico) em cada partido nas diversas legislaturas.
• Para isso vamos calcular as medianas e depois fazer um gráfico de linha.

• Calcular a mediana

median_dw1 <- congress %>%
  filter(party %in% c("Republican", "Democrat")) %>%
  group_by(party, congress) %>%
  summarize(median_dw1 = median(dwnom1)) %>%
  ungroup()

## `summarise()` has grouped output by 'party'. You can override using the
## `.groups` argument.

median_dw1

## # A tibble: 66 × 3
##    party    congress median_dw1
##    <chr>       <int>      <dbl>
##  1 Democrat       80     -0.126
##  2 Democrat       81     -0.207
##  3 Democrat       82     -0.179
##  4 Democrat       83     -0.174
##  5 Democrat       84     -0.223
##  6 Democrat       85     -0.231
##  7 Democrat       86     -0.259
##  8 Democrat       87     -0.258
##  9 Democrat       88     -0.285
## 10 Democrat       89     -0.292
## # ℹ 56 more rows

ggplot(data=median_dw1, aes(x=congress, y=median_dw1, color=party)) +
  geom_line() +
  labs(title = "DW-NOMINATE score (1st dimension)", x="Congress", y="Party") +
  theme_classic()

• O gráfico ilustra o distanciamento entre os dois partidos em termos de política econômica, os democratas ficam mais liberais (esquerda) e os republicanos mais conservadores (direita). A mudança dos republicanos é mais acentuada que a dos democratas.
• Esse fenômeno é chamado de polarização política.

• O que causou a polarização política? Essa é uma pergunta difícil de responder, o tema está agenda de pesquisa e gera debates na academia e outros meios. Alguns autores defendem que o aumento da desigualdade de renda está relacionado a essa polarização.
• Para medir a desigualdade de renda é comum usar o coeficiente de Gini.

• Para entender o coeficiente de Gini considere um gráfico onde o eixo horizontal representa o acumulado de pessoas em uma área específica (país, estado, etc) ordenado da menor para maior renda e o eixo vertical representa o acumulado de renda pertencente a indivíduos com renda igual ou menor que a da pessoa em um dado percentil de renda.
• A curva de Lorenz liga os pontos desse plano.
• Em uma distribuição de renda perfeitamente igualitária a curva de Lorenz é uma reta de 45o, isso ocorre porque x% da população detém, x% da renda da região. Essa reta é chamada de linha de igualdade.

• O coeficiente de Gini pode ser definido como a área entre a linha de igualdade e a curva de Lorenz dividida pela área abaixo da linha de igualdade. \[\small \mbox{coeficiente de Gini}=\frac{\mbox{área entre a linha de igualdade e a curva de Lorenz}}{\mbox{área abaixo da linha de igualdade}} \normalsize\]

• Na figura:

\[\mbox{coeficiente de Gini}=\frac{\mbox{área A}}{\mbox{área A + área B}}\]

• O coeficiente de Gini mede o grau de desigualdade de renda em uma sociedade. O coeficiente de Gini está entre zero (distribuição perfeitamente igualitária) e um (uma pessoa detém toda a renda).

• Para avaliar a relação entre desigualdade e polarização política comecemos comparando a evolução da medida de cada uma no tempo.
• Como medida de polarização usaremos a diferença entre as medianas dos dois partidos, como medida de desigualdade usaremos o coeficente de Gini.

polarization <- median_dw1 %>%
  pivot_wider(names_from = party,
              values_from = median_dw1) %>%
  mutate(polarization = Republican - Democrat)

polarization %>%
  ggplot(aes(congress, polarization)) +
  geom_point() +
  labs(title = "Politican polarization",
       x="Congress",
       y="Pepublican median -\n Democratic median") +
  theme_classic()

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

##   X year  gini
## 1 1 1947 0.376
## 2 2 1948 0.371
## 3 3 1949 0.378
## 4 4 1950 0.379
## 5 5 1951 0.363
## 6 6 1952 0.368

USGini %>%
  ggplot(aes(x=year, y=gini)) +
  geom_point() +
  labs(title = "Gini coefficent",
       x= "Year",
       y = "Income inequality") +
  theme_classic

• É fácil ver que as duas variáveis crescem no tempo, mas lembre isso pode acontecer por diversos motivos. Nunca esqueça: correlação não é causalidade!
• Por falar em correlação…

• Correlação, ou coeficiente de correlação, é uma das medidas mais usadas para representar relações entre duas variáveis. A correlação representa como, em média, duas variáveis se movimentam conjuntamente em relação às próprias médias.
• Grosso modo, a correlação diz se quando uma variável está acima (abiaxo) da média a outra também está.

• Antes de uma definição formal de correlação é útil definir o valor-z (z-score). O valor-z representa o número de desvios-padrão uma observação está acima ou abaixo da média. \[\mbox{valor-z de x_i}=\frac{x_i-\mbox{média de x}}{\mbox{desvio padrão de x}}\]

• O valor-z padroniza a variável de forma que tranformações lineares não mudam o valor-z de uma variável, por exemplo, mudanças na unidade de medida não mudam o valor-z. \[\mbox{valor-z de }(ax_i+b)=\frac{(ax_i+b)-\mbox{média de }(ax_i+b)}{\mbox{desvio-padrão de }(ax_i+b)}\] \[=\frac{a\times(x_i-{\mbox{média de }x})}{a\times\mbox{desvio-padrão de }x}\] \[=\mbox{valor-z de }x\]

• O valor-z da i-ésima observação de uma variável \(x\) mede o número de desvios padrões que a observação está abaixo ou acima da média. É definido por: \[\mbox{valor-z de }x_i=\frac{x_i-\bar{x}}{S_x}\] Onde \(\bar{x}\) e \(S_x\) são, respectivamente, a média e o desvio-padrão de \(x\). O valor-z é uma medida de desvio da média e não é modificado por transformações lineares em \(x\).

• Conhecido o valor-z, a correlação é definida por: \[\mbox{correlação}(x,y)=\frac{1}{n}\sum_{i=1}^n(\mbox{valor-z de }x_i \times \mbox{valor-z de }y_i)\]
• Assim como no desvio padrão, é comum usar \(n-1\) no lugar de \(n\) no denominador, continua valendo que essa mudança tem pouco efieto quando \(n\) é grande.

• Correlação, ou coeficiente de correlação, mede o grau em duas variáveis estão associadas uma a outra. É definido como: \[\mbox{correlação}(x,y)=\frac{1}{n}\sum_{i=1}^n\left(\frac{x_i-\bar{x}}{S_x}\times\frac{y_i-\bar{y}}{S_y}\right)\mbox{ ou }\\\frac{1}{n-1}\sum_{i=1}^n\left(\frac{x_i-\bar{x}}{S_x}\times\frac{y_i-\bar{y}}{S_y}\right)\] Onde \(\bar{x}\) e \(\bar{y}\) são as médias e \(S_x\) e \(S_y\) são os desvios-padrão das variáveis \(x\) e \(y\), respectivamente.Correlação vai de -1 a 1 e não depende da unidade de medida da variável.

• No R a função cor() calcula a correlação.
• Para calcular a correlação entre polarização política e coeficiente de Gini precisamos primeiro arrumar os dados, pois o coeficiente de Gini engloba todos os anos e cada legislatura dura por dois anos. Serão usados os coeficientes de Gini dos anos pares, correspondem ao segundo ano de cada legislatura.

gini_2yr <- USGini %>%
  filter(row_number() %% 2 ==0) %>%
  select(gini) %>%
  pull()

pol_annual <- polarization %>%
  select(polarization) %>%
  pull()

cor(gini_2yr, pol_annual)

## [1] 0.9418128

• A correlação é alta e positiva. Isso sugere que desigualdade de renda e polarização política andam juntas.
• Podemos falar de causalidade? Não, muitas variáveis mostram tendência de crescimento nesse período e, por óbvio, nem todas possuem relação de causalidade.

• Algumas vezes estamos interessados em comparar toda a distribuição de duas variáveis, não apenas medidas como média ou mediana, uma alternativa para fazer essa comparação é por meio do histograma de cada distribuição.
• Nesses tipos de comparação é importante que as escalas dos gráficos sejam as mesmas
• A função facet_grid() do ggplot2 permite fazer gráficos do tipo. No exemplo vamos comparar as distribuições da segunda dimensão do DW-NOMINATE score entre republicanos e democratas da 112o legislatura.

congress %>%
  filter(congress == 112, party %in% c("Republican", "Democrat")) %>%
  ggplot(aes(x=dwnom2, y=after_stat(density))) +
           geom_histogram(binwidth = .2) +
           facet_grid(party ~ .) +
           labs(x="Racial liberalism/conservatism dimension",
                y="Density") +
           theme_classic()

• As distribuições são semelhantes, mas é possível notar que a distribuição dos democratas tem uma cauda superior mais longa do que a dos republicanos (mais observações com valores mais altos).
• Outra diferença é que a distribuição dos republicanos está mais concentrada no centro do que a dos democratas.

• 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")

• jÁ vimos que a correlação representa o grau de associação entre duas variáveis. Uma correlação positiva signIfica 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).
• 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.

• 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>

• Curioso… 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 é de 97%!

• 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 e análise 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 lembre 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 tratamento, 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()