• Durante o curso vimos muitas formas de obter informações dos dados, mas como saber o quão certas são essas informações. Por exemplo, se a média do grupo de tratamento for diferente da média do grupo de controle, como saber se a diferença é por conta da amostra (lembram da variabilidade amostral?) ou se é por conta do tratamento? Em alguns exemplos dissemos que a diferença era grande, mas como definir esse grande? Enfim, como saber se o resultado observado vem do acaso?
• Responder, ainda que parcialmente, essas questões é o objetivo desta unidade. Para isso usaremos os conceitos de probabilidade que discutimos na unidade anterior.

• Por concistência com o livro esta unidade ficou com o nome de incerteza, no sentido que usaremos aqui incerteza faz referência a não termos certeza sobre o que é uma propriedade da amostra e o que é uma propriedade da população. Um problema relacionado à variabilidade amostral que tratamos na unidade anterior.
• Em outros contextos o termo incerteza é usado para descrever uma situação onde não conhecemos as distribuições, em contraponto com risco que ocorre quando conhecemos a distribuição. Usada nesse sentido é comum nos referimos a incerteza como incerteza knightiana.

• Qual o efeito de um aumento do gasto público na taxa de juros? Qual o efeito da redução da maioridade penal na taxa de assaltos? Qual o efeito de um aumento do bolsa família na oferta de trabalho? Qual o efeito de um aumento do salário mínimo no desemprego? Qual o efeito de barreiras à entrada na adoção de tecnologias? Qual o efeito da taxa de investimento no crescimento econômico?
• Perguntas como essas fazem parte da rotina de economistas aplicados.

• Em todos esses os casos usamos dados observados para estimar valores desconhecidos de quantidades que estamos interessados.
• Essas quantidades desconhecidas são chamadas de parâmetros, o método usado para estimar o parâmetro é chamado de estimador.

• A proporção de votos declarados em um candidato por quem respondeu uma pesquisa é um estimador para a proproção de votos do candidato entre o total de eleitores. A diferença entre a média do grupo de tratamento e a média do grupo de controle é um estimador para o efeito causal do tratamento.

• Quão bom é um estimador? É difícil responder essa pergunta, afinal não conhecemos o parâmetro que desejamos estimar.
• Porém, podemos usar de teoria ou simulação para avaliar o desempenho de um estimador em determinadas condições. Boa parte da sequência de econometria é dedicada a esse trabalho.

• Suponha que façamos uma pesquisa para estimar quem vai ganhar o segundo turno de uma eleição (ou afinal de um campeonato). Seja \(p\) a proporção de eleitores que vão votar no candidato de interesse, para seguir o livro suponha que o candidato de interesse é o Obama.

• Considere que a pesquisa foi respondida por \(n\) eleitores que foram escolhidos de forma independente. Cada eleitor é uma observação e será denotado por \(X_i\), se \(X_i =1\) o eleitor vota no candidato, se \(X_i=0\) o eleitor não vota no candidato. Supondo que existem apenas duas opções, \(\{X_1, X_2, \dots, X_n\}\) é uma sequência de variáveis aleatórias Bernoulli independentes e identicamente distribuídas.
• Nosso estimador para o parâmetro \(p\) será a média dos \(X_i's\), especificamente \(\bar{X}_n = \sum_{i=1}^n \frac{X_i}{n}\),

• O erro de estimação que pode ser uma medida que quão bom é um estimador é dado por: \[\mbox{erro de estimação} = \mbox{valor estimado}-\mbox{valor verdadeiro}=\bar{X}_n - p\]
• O problema com o erro estimação é que ele não pode ser calculado pois não conhecemos \(p\). Entretanto, podemos conhecer o valor esperado do erro de estimação.

• Considere um cenário hipotético onde fazemos a mesma pesquisa eleitoral infinitas vezes, cada pesquisa obteria um conjunto diferente de eleitores, todos com tamanho \(n\), e daria um diferente valor estimado para \(p\).
• Desta forma o erro de estimação seria diferente de uma amostra para outra, o erro de estimação é uma variável aleatória.

• A distribuição do erro amostral é chamada de distribuição amostral do estimador.
• No nosso exemplo o estimador é dado por \(\bar{X}_n = \sum_{i=1}^n \frac{X_i}{n}\), o númerador é uma soma de \(n\) Bernoulli i.i.d com probabilidade de sucesso igual a \(p\), portanto é uma binomial \((n,p)\), e o denominador é uma constante.
• Logo, \[E(\bar{X}_n)=\frac{1}{n}E(\sum_{i=1}^n X_i)=\frac{1}{n}np=p\]

• O valor esperado do erro de estimação é chamado de viés (bias), no nosso exemplo o viés é dado por: \[\mbox{viés}=E(\mbox{erro de estimação})\\=E(\bar{X}_n - p)=E(\bar{X}_n)-p=p-p=0\]
• Quando o viés é zero, como no caso acima, dizemos que o estimador é não-viesado.

• O resultado de que a média amostral é um estimador não-viesado da esperança não é exclusivo da distribuição Bernoulli. De fato, desde que a amostra seja construída por retiradas aleatórias e independetes da população, o resultado vale qualquer que seja a distribuição da população.

• Um estimador é dito não viesado se o valor esperado do estimador é igual ao parâmetro. Um estimador é dito consistente se converge para o parâmetro quando o tamanho da amostra cresce.
• Por exemplo, a média amostral, \(\bar{X}_n = \frac{\sum_{i=1}^nX_i}{n}\), é um estimador não viesado e consistente da média populacional (ou valor esperado), \(E(X)\), quando a amostra é obtida de forma aleatória e independente. \[E(\bar{X}_n)= E(X) \qquad \mbox{ e } \qquad \bar{X}_n \rightarrow E(X)\]

• O estimador de diferenças em médias usado para analisar controles aleatórios é não viesado para o efeito médio do tratamento, vamos ver a razão.
• Suponha que temos uma amostra de tamanho \(n\) na qual aplicaremos um experimento aleatorizado. O experimento consiste em um único tratamento, \(T_i\), que é igaual a um se a unidade \(i\) recebeu o tratamento e zero se a unidade \(i\) não recebeu o tratamento (grupo de controle).

• De forma aleatória escolhemos \(n_1\) unidades para receber o tratamento, as outras \(n-n_1\) unidades formam o grupo de controle.
• Essa forma de distribuir o tratamento, fixando com antecedência o número de unidades que receberão o tratamento, é chamada de aleatorização completa.
• Na aleatorização simples é feito um sorteio para cada unidade receber (ou não) o tratamento, os sorteios são independentes, por isso não sabemos com antecedência quantas unidades serão tratadas.

• Com aleatorização completa, existem \(_nC_{n_1}\) maneiras de distribuir \(n_1\) unidades para o grupo de tratamento e o restante para o grupo de controle. Cada uma dessas maneiras tem a mesma probabilidade, mas apenas uma acontece de fato.

• O primeiro parâmetro a considerar é o efeito de tratamento médio da amostra, SATE, que é definido por: \[\mbox{SATE}=\frac{1}{n}\sum_{i=1}^n \{ Y_i(1)-Y_i(0)\}\]
• Como vimos na segunda unidade, \(Y_i(1)\) representa o resultado se a unidade \(i\) receber o tratamento e \(Y_i(0)\) representa o resultado se a unidade \(i\) não recbeer o tratamento (estiver no grupo de controle). Não é possível observar \(Y_i(1)\) e \(Y_i(0)\), apenas um ou outro, desta forma não conhecemos o SATE.

• O estimador de diferenças em médias, \(\widehat{SATE}\), estima o SATE pela diferença entre a média dos tratados, \(T_i = 1\) e não tratados, \(T_i=0\). \[\widehat{SATE}= \mbox{média do grupo de tratamento}-\mbox{média do grupo de controle}\\=\frac{1}{n_1}\sum_{i=1}^n T_iY_i - \frac{1}{n-n_1}\sum_{i=1}^n (1-T_i)Y_i\]

• Para mostrar que \(E(\widehat{SATE})=SATE\), vamos usar que \(T_i\) é Bernoulli com valor esperado dado \(P(T_i=1)=\frac{n_1}{n}\)

\[E(\widehat{SATE})= E\left(\frac{1}{n_1}\sum_{i=1}^n T_iY_i(1)-\frac{1}{n-n_1}\sum_{i=1}^n (1-T_i)Y_i(0) \right)\\=\frac{1}{n_1} \sum_{i=1}^nE(T_i)Y_i(1)-\frac{1}{n-n_1}\sum_{i=1}^n E(1-T_i)Y_i(0)\\=\frac{1}{n_1}\sum_{i=1}^n \frac{n_1}{n}Y_i(1)-\frac{1}{n-n_1}\sum_{i=1}^n \left(1-\frac{n_1}{n}\right)Y_i(0)\\=\frac{1}{n_1} \sum_{i=1}^n\frac{n_1}{n}Y_i(1)-\frac{1}{n-n_1}\sum_{i=1}^n \left(\frac{n-n_1}{n}\right)Y_i(0)\\=\frac{1}{n}\sum_{i=1}^n\{Y_i(1)-Y_i(0) \}= \mbox{SATE}\]

• Suponha que você primeiro escolha aleatoriamente uma amostra de tamanho \(n\) de uma população e depois distribua aleatoriamente o tratamento entre as unidades da amostra como fizemos anteriormente. Com esses dois procedimentos podemos generalizar o resultado do experimento para a população.

• Para comprovar essa afirmação considere o efeito médio do tratamento na população, PATE (do inglês: population average treatment effect), dado por: \[\mbox{PATE}=E(Y_i(1)-Y_i(0))\]
• Por ter sido obtida de forma aleatória, a amostra é representativa da população. Desta forma, mesmo sendo o SATE não observado, o valor esperado do SATE é igual ao PATE. Portanto, como o estimador de diferenças em médias é não viesado para o SATE também é não viesado para o PATE.

• Em exprimentos aleatórios controlados, a diferença média entre os resultados dos grupos de tratamento e controle é um estimador não viesado do efeito médio do tratamento na amostra. Essa diferença também um estimador não viesado e consistente do efeito médio do tratamento na população, o resultado de consistência segue da lei dos grandes números aplicada aos grupos de tratamento e controle separadamente.

• Vamos usar uma simulação de Monte Carlo para ilustrar o resultado acima. Suponha que \(Y_i(0)\) é retirado de uma distribuição \(N(0,1)\) e que \(Y_i(1)\) é retirado de uma \(N(1,1)\), ou seja, \(Y(0) \sim N(0,1)\) e \(Y(1) \sim N(1,1)\).
• Dessa população é retirada aleatoriamente uma amostra de tamanho \(n\) e, de forma também aleatória, metade das unidades da amostra recebe o tratamento.
• Para cada unidade o efeito do tratamento é definido como \(\tau_i = Y_i(1)-Y_i(0)\)

• Repare que, como conhecemos os parâmetros da população, o PATE pode ser calculado: \[E(\tau_i)=E(Y_i(1)) - E(Y_i(0))=1-0-1\]
• O valor do SATE depende do sorteio das unidades que serão tratadas.

• Vamos simular esse experimento no R.

library(tidyverse)

n <- 100  
mu0 <- 0
sd0 <- 1
mu1 <- 1
sd1 <- 1

smpl <- tibble(id = seq_len(n),
               Y0 = rnorm(n, mean=mu0, sd=sd0),
               Y1 = rnorm(n, mean=mu1, sd=sd1),
               tau = Y1 - Y0)

SATE <- smpl %>% select(tau) %>%
  summarize(SATE = mean(tau)) %>% pull()

SATE

## [1] 0.7475618

• Podemos repetir o exercício acima várias vezes mudando as unidades que recebem o tratamento, para isso convém escrever uma função.

sim_treat <- function(smpl) {
  n <- nrow(smpl)
  idx <- sample(seq_len(n), floor(nrow(smpl)/2), replace = FALSE)
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>% 
    pivot_wider(names_from = treat, values_from = Y_obs) %>%
    rename(Y1_mean = `1`, Y0_mean = `0`) %>%
    mutate(diff_mean = Y1_mean - Y0_mean,
           est_error = diff_mean - SATE)
}

sim_treat(smpl)

## # A tibble: 1 × 4
##   Y0_mean Y1_mean diff_mean est_error
##     <dbl>   <dbl>     <dbl>     <dbl>
## 1   0.148   0.723     0.575    -0.172

• Antes de seguir com exemplo, é válido parar para entender melhor a função que acabamos de criar. Para isso vamos “quebrar” a função e ver o que faz cada pedaço.

n <- nrow(smpl)
idx <- sample(seq_len(n), floor(nrow(smpl)/2), replace = FALSE)

n

## [1] 100

idx

##  [1]  24  70   1  13  65  17  52  42  18  84  37  38  19   9  41  33  59  74  21
## [20]  87  79  93  90  91  28  80  67  76  71  35  45  62  69  46  95  20  50  99
## [39]  51  66  86  88  56  58  75  26   6 100  97  23

smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)

head(smpl)

## # A tibble: 6 × 5
##      id     Y0     Y1    tau treat
##   <int>  <dbl>  <dbl>  <dbl> <int>
## 1     1  1.27   1.93   0.654     1
## 2     2 -1.18   0.864  2.04      0
## 3     3  0.386  0.612  0.226     0
## 4     4  2.31  -0.936 -3.25      0
## 5     5 -0.274  1.54   1.81      0
## 6     6 -1.29   1.17   2.46      1

smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0))

## # A tibble: 100 × 6
##       id     Y0     Y1    tau treat  Y_obs
##    <int>  <dbl>  <dbl>  <dbl> <int>  <dbl>
##  1     1  1.27   1.93   0.654     1  1.93 
##  2     2 -1.18   0.864  2.04      0 -1.18 
##  3     3  0.386  0.612  0.226     0  0.386
##  4     4  2.31  -0.936 -3.25      0  2.31 
##  5     5 -0.274  1.54   1.81      0 -0.274
##  6     6 -1.29   1.17   2.46      1  1.17 
##  7     7 -0.463  0.454  0.917     0 -0.463
##  8     8  1.07   1.24   0.161     0  1.07 
##  9     9  0.537  0.430 -0.107     1  0.430
## 10    10  0.137 -1.09  -1.23      0  0.137
## # ℹ 90 more rows

smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs))

## # A tibble: 2 × 2
##   treat  Y_obs
##   <int>  <dbl>
## 1     0 0.0679
## 2     1 0.749

smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>% 
    pivot_wider(names_from = treat, values_from = Y_obs)

## # A tibble: 1 × 2
##      `0`   `1`
##    <dbl> <dbl>
## 1 0.0679 0.749

smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>% 
    pivot_wider(names_from = treat, values_from = Y_obs) %>%
    rename(Y1_mean = `1`, Y0_mean = `0`) %>%
    mutate(diff_mean = Y1_mean - Y0_mean,
           est_error = diff_mean - SATE)

## # A tibble: 1 × 4
##   Y0_mean Y1_mean diff_mean est_error
##     <dbl>   <dbl>     <dbl>     <dbl>
## 1  0.0679   0.749     0.681   -0.0669

• Agora que entendemos a função que criamos vamos seguir com a simulação. Usaremos afunção map_df() para aplicar a função 500 vezes na nossa amostra e observar o erro de estimação.

sims <- 500
sate_sims <- map_df(seq_len(sims), ~ sim_treat(smpl))

summary(sate_sims$est_error)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.3789586 -0.0927292  0.0082436  0.0006646  0.0825342  0.4224462

• Podemos usar uma estratégia semelhante para analisar o estimador em relação ao PATE.
• Para isso a função deve ser modificiada de forma a “sortear” uma amostra da população antes de aplicar os procedimentos que estão na função sim_treat().

• Lembre que o PATE é dado por:

PATE <- mu1 - mu0

sim_pate <- function(n, mu0, mu1, sd0, sd1){
  PATE = mu1 - mu0
  smpl <- tibble(Y0 = rnorm(n ,mean=mu0, sd=sd0),
                 Y1 = rnorm(n, mean=mu1, sd=sd1),
                 tau = Y1 - Y0)
  idx <- sample(seq_len(n), floor(nrow(smpl)/2), replace = FALSE)
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>% 
    pivot_wider(names_from = treat, values_from = Y_obs) %>%
    rename(Y1_mean = `1`, Y0_mean = `0`) %>%
    mutate(diff_mean = Y1_mean - Y0_mean,
           est_error = diff_mean - PATE)
}

sims <- 500
pate_sims <- map_df(seq_len(sims), ~ sim_pate(n, mu0, mu1, sd0, sd1))
summary(pate_sims$est_error)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.549557 -0.118001  0.002255  0.003158  0.135687  0.466888

• Nos dois casos as médias dos erros de estimação ficaram próximas de zero, o erro de estimação em relação ao PATE tem mais variabilidade porque são duas etapas aleatórias: seleção da amostra e seleção das unidades tratadas, isso coloca mais ruído.

• Um estimador não viesado, mas com muita variação pode ser pouco útil na prática. Tal estimador acertaria o parâmetro na média, mas as realizações, que é o que temos, poderiam estar muito longe do parâmetro.
• O estimador de diferenças das médias é não viesado, mas os valores estimados mudam muito de uma simulação para outra.

pate_sims %>%
  ggplot(aes(diff_mean, after_stat(density))) +
  geom_histogram(binwidth = 0.1) +
  geom_vline(xintercept = PATE, color="black", linewidth=0.5) +
  labs(title = "Sampling distribution",
       x="Difference-in-means estimator",
       y="Density") +
  theme_classic()

• Uma maneira de medir a variação do nosso estimador é calcular o desvio padrão do valor obtido para o estimador nas diversas amostra criadas.
• O desvio padrão mostra a distância média do estimador para a média do próprio estimador.

pate_sims %>%
  select(diff_mean) %>%
  summarize(sd = sd(diff_mean))

## # A tibble: 1 × 1
##      sd
##   <dbl>
## 1 0.184

• Como o estimador de diferenças das médias é não viesado, a média do estimador é igual ao SATE. Desta forma, o desvio padrão da média das amostras é igual a raiz do erro quadrado médio (RMSE).

RMSE <- pate_sims %>%
  summarize(rmse = sqrt(mean(est_error^2)))

• Se o estimador for viesado, o desvio padrão da distribuição amostral do estimador pode ser diferente do RMSE. Especificamente, se o estimador for viesado, o erro médio quadrado (MSE), corresponde ao quadrado do RMSE, é igual à variância do estimador mais o quadrado do viés.

\[\mbox{MSE}=E[(\hat{\theta}-\theta)^2]=E[\hat{\theta}-E(\hat{\theta})+E(\hat{\theta})-\theta]^2\\=E[(\hat{\theta}-E(\hat{\theta}))+(E(\hat{\theta})-\theta)]^2\\=E[(\hat{\theta}-E(\hat{\theta}))^2 +2E(\hat{\theta}-E(\hat{\theta}))(E(\hat{\theta})-\theta)+(E(\hat{\theta})-\theta)^2]\\=E[(\hat{\theta}-E(\hat{\theta}))] + E[(E(\hat{\theta})-\theta)^2]=\mbox{variância} + \mbox{viés}^2\]

• Para chegar ao resultado usamos que \(2E(\hat{\theta}-E(\hat{\theta}))(E(\hat{\theta})-\theta)=0\). Para verificar que isso é verdade note que:

\[E(\hat{\theta}-E(\hat{\theta}))(E(\hat{\theta})-\theta)=E[\hat{\theta}E(\hat{\theta})-\hat{\theta}\theta-E(\hat{\theta})^2+E(\hat{\theta})\theta]\\=E(\hat{\theta})^2-E(\hat{\theta})\theta-E(\hat{\theta})^2+E(\hat{\theta})\theta=0\]

• Uma lição importante da fórmula \(\mbox{MSE}=\mbox{variância} + \mbox{viés}^2\) é que um estimador não viesado pode ter uma MSE maior do que a de um estimador viesado a depender da variãncia de cada um.

• Até agora medimos variação de um estimador usando o desvio padrão da distribuição amostral, o problema é que essa medida não pode ser obtida diretamente dos dados. Em aplicações só temos uma amostra do processo gerador de dados. Por isso recorremos ao erro padrão que é uma estimativa do desvio padrão da distribuição amostral do estimador.
• O erro padrão pode ser obtido diretamente dos dados.

• Para mensurar a variabilidade de um estimador podemos usar o erro padrão, que é uma estimativa do desvio padrão da distribuição amostral. Uma medida de acurácia é a raiz do erro quadrado médio (RMSE), que mede o desvio médio de um estimador, \(\hat{\theta}\), em relação ao parâmetro, \(\theta\). O quadrado do RMSE, chamado de erro quadrado médio (MSE) é igual a soma da variância com o quadrado do viés do estimador. \[E(\hat{\theta}-\theta)^2=V(\hat{\theta})+[E(\hat{\theta}-\theta)]^2\]

• No exemplo de estimar a proproção de votos para Obama tínhamos uma única amostra aleatória de tamanho \(n\). Podemos interpretar a resposta de cada entrevistado como uma realização i.i.d. de uma Bernoulli com probabilidade de sucesso igual a \(p\).

• Usamos a média, que representa a proporção de sucessos, como estimador, \(\bar{X}_n=\frac{\sum_{i=1}^n X_i}{n}\), a variância do estimador é dada por: \[V(\bar{X}_n)=\frac{1}{n^2}V(\sum_{i=1}^n X_i)=\frac{1}{n^2}\sum_{i=1}^nV(X_i)\\=\frac{np(1-p)}{n^2}=\frac{p(1-p)}{n}\]

• Para um dado valor de \(n\), a variância do estimador é função de \(p\)

• Não conhecemos \(p\), mas sabemos que \(\bar{X}_n\) é um estimador não viesado e consistente de \(p\), que pode ser usado para estimar o erro padrão como a raiz da variância.

\[\mbox{erro padrão}=\sqrt{\frac{\bar{X}_n(1-\bar{X}_n)}{n}}\]

• Suponha que temos uma amostra de \(n\) variáveis aleatórias, \(\{X_1,X_2,\dots,X_n\}\), i.i.d. O erro padrão da média amostral \(\bar{X}_n=\frac{\sum_{i=1}^n X_i}{n}\) é dado por:

\[\mbox{erro padrão da média amostral}=\sqrt{\widehat{V(\bar{X}_n)}}=\sqrt{\frac{\widehat{V(X)}}{n}}\]

• Para estimar a variância da população \(V(X)\) usamos a variância da amostra \(\frac{\sum_{i=1}^n(X_i - \bar{X}_n)^2}{n-1}\). O númerador é \(n-1\) no lugar de \(n\) porque para estimar a variância precisamos estimar a média o que resulta na perda de um grau de liberdade.

• Suponha que temos uma amostra i.i.d. de tamanho \(n\), \(\{X_1,X_2,\dots,X_n\}\), e uma amostra i.i.d. de tamanho \(m\), \(\{Y_1,Y_2,\dots,Y_n\}\). O erro padrão do estimador de diferenças da média, \(\frac{\sum_{i=1}^n X_i}{n}-\frac{\sum_{i=1}^m Y_i}{m}\), é dado por: \[\mbox{erro padrão da diferença das médias}=\sqrt{\frac{\widehat{V(X)}}{n}+\frac{\widehat{V(Y)}}{m}}\]

• Vamos usar de Monte Carlo para checar esse último resultado.

sim_pate_se <- function(n, mu0, mu1, sd0, sd1){
  PATE = mu1 - mu0
  smpl <- tibble(Y0 = rnorm(n ,mean=mu0, sd=sd0),
                 Y1 = rnorm(n, mean=mu1, sd=sd1),
                 tau = Y1 - Y0)
  idx <- sample(seq_len(n), floor(nrow(smpl)/2), replace = FALSE)
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  smpl %>%
    mutate(Y_obs = if_else(treat==1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(mean = mean(Y_obs),
              var = var(Y_obs),
              nobs=n()) %>%
    summarize(diff_mean = diff(mean),
              se=sqrt(sum(var/nobs)),
              std_error = diff_mean-PATE)
}

• Uma simulação para testar a função

sim_pate_se(n, mu0, mu1, sd0, sd1)

## # A tibble: 1 × 3
##   diff_mean    se std_error
##       <dbl> <dbl>     <dbl>
## 1      1.06 0.191    0.0598

• 500 simulações

sims <- 500
pate_sims_se <- map_df(seq_len(sims), 
                       ~ sim_pate_se(n, mu0, mu1, sd0, sd1))

sd_se <- pate_sims_se %>%
  summarize(sd = sd(diff_mean),
            mean_se = mean(se))

sd_se

## # A tibble: 1 × 2
##      sd mean_se
##   <dbl>   <dbl>
## 1 0.201   0.199

• Como conhecemos os processos geradores de dados de \(Y(1)\) e \(Y(0)\), podemos calcular o desvio padrão do erro analiticamente:

\[\sqrt{\frac{V(Y(1))}{n_1}+\frac{V(Y(2))}{n-n_1}}=\sqrt{\frac{1}{50}+\frac{1}{50}}\\=\sqrt{\frac{1}{25}}=\frac{1}{5}\]

• Já vimos como encontrar a média e a variância da dsitribuição amostral de um estimador, agora queremos saber mais sobre a distribuição propriamente dita.

• Se o processor gerador de dados for \(N(\mu, \sigma^2)\), não é difícil encontrar a distribuição de \(\bar{X}_n=\frac{\sum_{i=1}^n X_i}{n}\). Basta lembrar que se \(\{X_1, X_2, \dots, X_n\}\) são i.i.d. \(N(\mu,\sigma^2)\), então: \[\sum_{i=1}^n X_i \sim N(n \mu, n\sigma^2)\]
• Logo, \[\frac{\sum_{i=1}^n X_i}{n} \sim N\left(\mu, \frac{\sigma^2}{n}\right)\]

• O problema é que nem sempre o processo gerador de dados segue uma distribuição normal, fica pior, via de regra nem sabemos qual é a distribução do processo gerador de dados.
• Mesmo sem conhecer a distribuição do processo gerador de dados, podemos determinar a distribuição de \(\bar{X}_n=\frac{\sum_{i=1}^n X_i}{n}\).

• Para isso lembre que o teorema central do limite diz que se \(\{X_1, X_2, \dots, X_n\}\) são i.i.d. com média \(E(X)\) e variância \(V(X)\), então: \[\frac{\bar{X}_n- E(X)}{\sqrt{\frac{V(X)}{n}}} \quad \mbox{converge em distribuição para} \quad N(0,1)\]

• Desta forma, multiplicando \(\frac{\bar{X}_n- E(X)}{\sqrt{\frac{V(X)}{n}}}\) por \(\sqrt{\frac{V(X)}{n}}\):

\[\bar{X}_n- E(X) \sim N\left(0, \frac{V(x)}{n}\right)\]

• Somando \(E(X)\):

\[\bar{X}_n \sim N\left(E(X), \frac{V(x)}{n}\right)\]

• Por exemplo, se \(\{X_1, X_2, \dots, X_n\}\) são variáveis aleatórias de Bernoulli com probabilidade de sucesso \(p\), então, para um \(n\) suficientemente grande, \(\bar{X}_n \sim N\left(p, \frac{p(1-p)}{n}\right)\).

• Esse resultado nos permite construir uma medida de incerteza chamada intervalo de confiança, que está relacionado à probabilidade do valor do parâmetro estar em um determinado intervalo.

• Para construir um intervalo de confiança primeiro é preciso definir o nível de confiança, que costuma ser visto como a probabilidade do parâmetro estar no intervalo, embora, estritamente falando, represente a percentagem de intervalos gerados em repetidas simulações que contém o verdadeiro valor do parâmetro.
• Por exemplo, se o nível de confiança for de 95% isso quer dizer que 95% dos intervalos gerados em “infitas” simulações contém o verdadeiro valor do parâmetro.

• É comum definir o nível de confiança como \((1-\alpha)\times 100\%\) onde \(\alpha\) é um número entre zero e um. Desta forma, um nível de confiança de 95% é exepresso por \(\alpha=0,\!05\).

• Formalmente o intervalo de confiança para média com nível de confiança dado por \(\alpha\) é dado por: \[\mbox{CI}(\alpha)=[\bar{X}_n - z_{\frac{\alpha}{2}}\times \mbox{erro padrão}, \bar{X}_n + z_{\frac{\alpha}{2}}\times \mbox{erro padrão}]\]
• Onde \(z_{\frac{\alpha}{2}}\) é chamado de valor crítico e corresponde ao quantil \((1-\frac{\alpha}{2})\) da distribuição normal padrão, ou seja, \(P(Z>\frac{\alpha}{2})=1-P(Z\leq \frac{\alpha}{2})=1-\frac{\alpha}{2}\), onde \(Z \sim N(0,1)\).

• A próxima figura ilustra a relação entre o valor crítico, \(z\), e o nível de confiança.
• A área total é um, as áreas cinzas são dadas por \(\frac{\alpha}{2}\) cada, de forma que a área verde é \(1-\alpha\).
• Desta forma, se desejamos um nível de confiança de 95%, temos de escolher um \(\alpha\) tal que as áreas cinzas sejam 2,5% cada.

• A função qnorm() permite calcular o valor crítico para cada nível de confiança.
• Essa função calcula a inversa da CDF de uma distribuição normal padrão, o argumento da função tem uma probbailidade \(p\) e a função calcula o quantil tal que \(p=P(X \leq q)\).

#alpha = 0,01 -> nível de confiança 99%
qnorm(0.995) # 1-0.01/2

## [1] 2.575829

#alpha = 0,05 -> nível de confiança 95%
qnorm(0.975) # 1-0.05/2 

## [1] 1.959964

#alpha = 0,1 -> nível de confiança 90%
qnorm(0.95) # 1-0.1/2 

## [1] 1.644854

• No exemplo da pesquisa eleitoral, se 600 de 1000 eleitores declaram voto em Obama então \(\bar{X}_n=0,\!6\) e o erro padrão é \(\sqrt{\frac{0,\!6 \times (1-0,\!6)}{1000}}\).

n <- 1000
x_bar <- 0.6

se <- sqrt(x_bar * (1-x_bar)/n)
levels <- c(0.99, 0.95, 0.9)

tibble(level = levels) %>%
  mutate(ci_lower = x_bar - qnorm(1-(1-level)/2)*se,
         ci_upper = x_bar + qnorm(1-(1-level)/2)*se)

## # A tibble: 3 × 3
##   level ci_lower ci_upper
##   <dbl>    <dbl>    <dbl>
## 1  0.99    0.560    0.640
## 2  0.95    0.570    0.630
## 3  0.9     0.575    0.625

• O intervalo de confiança de um estimador \(\hat{\theta}\) pode ser obtido usando os seguintes passos:
  • Escolha o nível de confiança \((1-\alpha)\times 100\%\) especificando o valor de \(\alpha\), a esolha mais comum é \(\alpha=0,\!05\) que equivale a um nível de confiaça de 95%.
  • Calcule a média e a variãncia amostral do estimador.
  • Calcule o erro padrão.
  • Calcule o valor crítico.
  • Calcule o limite inferior e o limite superior.

• Vamos usar as funções criadas para simular o PATE para ilustrar o signficado de intervalo de confiança.

level <- 0.95

pate_sims_ci <- pate_sims_se %>%
  mutate(ci_lower = diff_mean - qnorm(1 - (1-level)/2)*se,
         ci_upper = diff_mean + qnorm(1 - (1-level)/2)*se,
         includes_pate = PATE > ci_lower & PATE < ci_upper)

glimpse(pate_sims_ci)

## Rows: 500
## Columns: 6
## $ diff_mean     <dbl> 0.8218692, 1.0451775, 1.0508220, 1.0515216, 1.3311669, 1…
## $ se            <dbl> 0.1923786, 0.2052930, 0.2013800, 0.2121763, 0.1819019, 0…
## $ std_error     <dbl> -0.17813084, 0.04517751, 0.05082202, 0.05152156, 0.33116…
## $ ci_lower      <dbl> 0.4448140, 0.6428107, 0.6561245, 0.6356637, 0.9746458, 0…
## $ ci_upper      <dbl> 1.1989243, 1.4475443, 1.4455195, 1.4673794, 1.6876880, 1…
## $ includes_pate <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, …

• Qual o percentual dos intervalos que incluem o valor correto do PATE?

pate_sims_ci %>%
  summarize(coverage = mean(includes_pate))

## # A tibble: 1 × 1
##   coverage
##      <dbl>
## 1     0.95

• Criar uma função para repetir o exercício acima para diferentes níveis de confiança.

pate_sim_coverage <- function(.data, level = 0.95){
  mutate(.data,
         ci_lower = diff_mean - qnorm(1 - (1-level)/2)*se,
         ci_upper = diff_mean + qnorm(1 - (1-level)/2)*se,
         includes_pate = PATE > ci_lower & PATE < ci_upper) %>%
    summarize(coverage = mean(includes_pate))
}

• Testar a função, nível de confiança de 95%

pate_sim_coverage(pate_sims_se)

## # A tibble: 1 × 1
##   coverage
##      <dbl>
## 1     0.95

pate_sim_coverage(pate_sims_se, level = 0.95)

## # A tibble: 1 × 1
##   coverage
##      <dbl>
## 1     0.95

• Testar a função, nível de confiança de 99%

pate_sim_coverage(pate_sims_se, level = 0.99)

## # A tibble: 1 × 1
##   coverage
##      <dbl>
## 1     0.99

• Testar a função, nível de confiança de 90%

pate_sim_coverage(pate_sims_se, level = 0.9)

## # A tibble: 1 × 1
##   coverage
##      <dbl>
## 1    0.896

• Para ilustrar a relação entre nível de confiança e tamanho da amostra vamos retomar o exemplo da pesquisa eleitoral.
• Vamos começar construindo uma função que, dados \(n\) e \(p\), calcula o intervalo de confiança e verifica se o verdadeiro valor do parâmetro, \(p\), está no intervalo.

binom_ci_contains <- function(n, p, alpha=0.05){
  x <- rbinom(n, size=1, prob = p)
  x_bar <- mean(x)
  se <- sqrt(x_bar * (1-x_bar)/n)
  ci_lower <- x_bar - qnorm(1 - alpha/2)*se
  ci_upper <- x_bar + qnorm(1 - alpha/2)*se
  (ci_lower <= p) & (p <= ci_upper)
}

p <- 0.6
n <- 10
binom_ci_contains(n=n, p=p, alpha=0.05)

## [1] TRUE

• A função map_lgl() é semelhante à função map_df() porém retorna um objeto do tipo logical em vez de um data.frame. Vamos usar essa função para aplicar várias vezes a função binom_ci_contais() com amostras de tamanho \(n\).

mean(map_lgl(seq_len(sims), ~ binom_ci_contains(n,p)))

## [1] 0.894

• Para repetir o exercício com diferentes valores de \(n\) vamos colocar o código acima em uma função.

binom_ci_coverage <- function(n, p, sims){
  mean(map_lgl(seq_len(sims), ~ binom_ci_contains(n,p)))
}

tibble(n=c(10, 100, 1000)) %>%
  mutate(coverage = map_dbl(n, binom_ci_coverage, p=p, sims=sims))

## # A tibble: 3 × 2
##       n coverage
##   <dbl>    <dbl>
## 1    10    0.876
## 2   100    0.946
## 3  1000    0.944

• Em pesquisas eleitorais é comum que os resultados venham acompanhados da margem de erro, que é calculada a partir do intervalo de confiança. Quando dizemos que um candidato tem 60% das intenções de voto com margem de erro de3% para mais ou para menos, em geral estamos dizendo que o intervalo de confiança de 95% é [57, 63]. \[\mbox{margem de erro}= \pm 1,\!96 \times \sqrt{\frac{\bar{X}_n(1-\bar{X}_n)}{n}}\]

• Para um dado \(n\), o maior valor da margem de erro ocorre quando \(p=0,\!5\). Esse valor é dado por: \[\pm 1,\!96 \times \sqrt{\frac{0,\!5 \times 0,\!5}{n}} = \pm 1,\!96 \times \sqrt{\frac{0,\!25}{n}}\\=\pm 1,\!96 \times 0,\!5\times \frac{1}{\sqrt{n}} \approx \pm \frac{1}{\sqrt{n}}\]

• Dessa relação vem uma famosa regra de bolso para calcular tamanho de amostras. A regra diz que calculando com a maior margem de erro, o tamanho da amostra é dado por \(n \approx \frac{1}{\mbox{margem de erro}^2}\).
• Por exemplo, para uma margem de erro de \(\pm 3\) devemos usar uma amostra de tamanho \(n=\frac{1}{0,\!03^2} \approx 1111\) observações.

• A margem de erro da proporção estimada em pesquisas, \(\bar{X}_n\) se refere ao meio tamanho do intervalo de confiança de 95%, ou \(z_{0,\!025}\times \mbox{erro padrão}=1,\!96 \times \sqrt{\frac{\bar{X}_n(1-\bar{X}_n)}{n}}\). A relação entre margem de erro e tamanho da amostra é dada por: \[n = \frac{1,\!96^2 p (1-p)}{\mbox{margem de erro}^2}\]
• A fórumla pode ser usada para determinar o tamanho da amostra necessária para pesquisa.

• Podemos calcular o tamanho da amostra para diferentes níveis de margem de erro, para ilsutrar essa relaçao vamos criar uma função relacionando \(n\) e \(p\) para diferentes margens de erro.

moe_pop_prop <- function(MoE) {
  tibble(p = seq(from = 0.01, to = 0.99, by=0.01),
         n = 1.96^2 * p *(1-p)/MoE^2,
         MoE = MoE)
}

glimpse(moe_pop_prop(0.01))

## Rows: 99
## Columns: 3
## $ p   <dbl> 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, …
## $ n   <dbl> 380.3184, 752.9536, 1117.9056, 1475.1744, 1824.7600, 2166.6624, 25…
## $ MoE <dbl> 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, …

MoE <- c(0.01, 0.03, 0.05)
props <- map_df(MoE, moe_pop_prop)

props %>%
  ggplot(aes(x=p, y=n, color = factor(MoE))) +
  geom_line(aes(linetype = factor(MoE))) +
  labs(color = "Margin of error",
       linetype = "Margin of error",
       x = "Population proportion",
       y = "Sample size") +
  theme_classic() +
  theme(legend.position = "bottom") +

• Tamanho da amosta com \(p=0,\!05\) e diferentes margens de erro.

props %>%
  filter(p == 0.5) %>%
  select(MoE,n)

## # A tibble: 3 × 2
##     MoE     n
##   <dbl> <dbl>
## 1  0.01 9604 
## 2  0.03 1067.
## 3  0.05  384.

• Na Unidade 4 comparamos os valores previstos e realizados para a proporção de votos em Obama para cada estado nas eleições de 2008.
• Vamos repetir o exercício, porém agora vamos considerar as margens de erro. Os dados estão em pres08, suponha que tamanho de cada amostra é 1000.

library(lubridate)
pres08 <- read.csv("pres08.csv")
polls08 <- read.csv("polls08.csv") %>%
  mutate(middate = as.Date(middate))

ELECTION_DATE <- ymd(20081104)

polls08 <- polls08 %>%
  mutate(DaysToElection = as.integer(ELECTION_DATE - middate))

poll_pred <- polls08 %>%
  group_by(state) %>%
  filter(DaysToElection == min(DaysToElection)) %>%
  summarize(Obama = mean(Obama)/100)

head(poll_pred)

## # A tibble: 6 × 2
##   state Obama
##   <chr> <dbl>
## 1 AK    0.39 
## 2 AL    0.36 
## 3 AR    0.44 
## 4 AZ    0.465
## 5 CA    0.6  
## 6 CO    0.52

• Acrescentar os intervalos de confiança

sample_size <- 1000
alpha <- 0.05

poll_pred <- poll_pred %>%
  mutate(se = sqrt(Obama * (1-Obama)/sample_size),
         ci_lwr = Obama + qnorm(alpha/2)*se,
         ci_upr = Obama + qnorm(1-alpha/2)*se)

head(poll_pred)

## # A tibble: 6 × 5
##   state Obama     se ci_lwr ci_upr
##   <chr> <dbl>  <dbl>  <dbl>  <dbl>
## 1 AK    0.39  0.0154  0.360  0.420
## 2 AL    0.36  0.0152  0.330  0.390
## 3 AR    0.44  0.0157  0.409  0.471
## 4 AZ    0.465 0.0158  0.434  0.496
## 5 CA    0.6   0.0155  0.570  0.630
## 6 CO    0.52  0.0158  0.489  0.551

• Adicionar a proporção de votos que Obama teve em cada estado e checar se está no intervalo.

poll_pred <- poll_pred %>% 
  left_join(select(pres08, state, actual = Obama), by="state") %>%
  mutate(actual = actual/100,
         covers = (ci_lwr <= actual) & (actual <= ci_upr))

• Ilustrar o resultado com um gráfico

poll_pred %>%
  ggplot(aes(actual, Obama, color = as.factor(covers))) +
  geom_abline(intercept = 0, slope = 1, color = "black", linewidth = 0.5) +
  geom_pointrange(aes(ymin=ci_lwr, ymax=ci_upr), shape=1) +
  scale_y_continuous(limits = c(0,1)) +
  scale_x_continuous(limits = c(0,1)) +
  scale_color_discrete("CI includes result?") +
  labs(x = "Obama'a vote share",
       y= "Poll prediction") +
  coord_fixed() +
  theme_classic()

• Proporção dos intervalos contendo os valores que de fato ocorreram

poll_pred %>%
  summarize(mean(covers))

## # A tibble: 1 × 1
##   `mean(covers)`
##            <dbl>
## 1          0.588

• O resultado ficou muito abaixo do esperado, apenas 58,8% continham o valor obervado para o parâmetro
• Uma explicação para isso é que a estimativa de votos para Obama estava viesada, se isso acontece os intervalos de confiança são deslocados para cima (ou para baixo) e o erro padrão é alterado.

• Para avaliar essa possibilidade vamos calcular o viés e corrigir os intervalos, repare que isso só pode ser feito porque conhecemos os resultados da eleição.

poll_pred <- poll_pred %>%
  mutate(bias = Obama - actual) %>%
  mutate(Obama_bc = Obama - mean(bias),
         se_bc = sqrt(Obama_bc * (1-Obama_bc)/ sample_size),
         ci_lwr_bc = Obama_bc + qnorm(alpha/2)*se_bc,
         ci_upr_bc = Obama_bc + qnorm(1-alpha/2)*se_bc,
         covers_bc = (ci_lwr_bc <= actual) & (actual <= ci_upr_bc))

poll_pred %>% summarize(mean(bias))

## # A tibble: 1 × 1
##   `mean(bias)`
##          <dbl>
## 1      -0.0268

mean(poll_pred$bias)

## [1] -0.02679739

• Proporção Proporção dos intervalos contendo os valores que de fato ocorreram com correção pelo viés

poll_pred %>%
  summarize(mean(covers_bc))

## # A tibble: 1 × 1
##   `mean(covers_bc)`
##               <dbl>
## 1             0.765

• Melhorou um bocado, embora ainda esteja bem abaixo dos 95%.

• O projeto STAR (Student-Teacher Achievment Ratio) foi conduzido na década de 1980 para avaliar o impacto do tamanho das salas no desempenho dos estudantes em testes.
• Para isso os estudantes foram distribuídos de forma aleatória em salas pequenas, salas de tamanho regular e salas de tamanho regular com ajudante.

• Um estudo sobre o projeto STAR pode ser encontrado em: Frederick Mosteller. The Tennessee study of class size in the early school grades. Future of Children, v.5, n.2, 1995.
• Os dados estão em STAR.

• Variáveis:
  • race: raça do estudante (branco=1, afro-americano-2, asiático=3, hispânico=4, nativo americano=5, outros=6)
  • classtype: tipo de classe (pequena=1, regular=2, regular com ajudante=3)
  • g4math: nota da parte de matemática do exame padrão para estudantes da fourth-grade
  • g4reading: nota da parte de leitura do exame padrão para estudantes da fourth-grade
  • yearssmall: número de anos em classes pequenas
  • hsgrad: se terminou a *high-school** (terminou=1, não terminou=0)

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

STAR <- STAR %>%
  mutate(classtype = factor(classtype, 
                            labels = c("Small class", "Regular class", 
                                       "Regular class with aide")))

glimpse(STAR)

## Rows: 6,325
## Columns: 6
## $ race       <int> 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1,…
## $ classtype  <fct> Regular class with aide, Regular class with aide, Regular c…
## $ yearssmall <int> 0, 0, 0, 4, 0, 0, 4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0,…
## $ hsgrad     <int> NA, NA, 1, NA, NA, NA, NA, NA, 1, 1, NA, NA, 1, NA, 1, NA, …
## $ g4math     <int> NA, 706, 711, 672, NA, NA, 668, NA, 709, 698, NA, NA, NA, N…
## $ g4reading  <int> NA, 661, 750, 659, NA, NA, 657, NA, 725, 692, NA, NA, NA, N…

• Média da nota em leitura por tipo de classe

classtype_means <- STAR %>%
  group_by(classtype) %>%
  summarize(g4reading = mean(g4reading, na.rm=TRUE))

classtype_means

## # A tibble: 3 × 2
##   classtype               g4reading
##   <fct>                       <dbl>
## 1 Small class                  723.
## 2 Regular class                720.
## 3 Regular class with aide      721.

• Agora vamos olhar as distribuições.

classtypes_used <- c("Small class", "Regular class")

STAR %>%
  filter(classtype %in% classtypes_used, !is.na(g4reading)) %>%
  ggplot(aes(g4reading, after_stat(density))) +
  geom_histogram(binwidth = 20) +
  geom_vline(data = filter(classtype_means, classtype %in% classtypes_used),
             mapping = aes(xintercept=g4reading),
             color = "blue", linewidth = 0.5) +
  facet_grid(classtype ~.) +
  labs(x="Fourth-grade reading test scores",
       y="Density") +
  theme_classic()

• As médias parecem próximas uma da outra e as distribuições apresentam padrões semelhantes. Será que o tamanho da classe faz diferença?
• Para tentar responder essa pergunta é interessante construir os intervalos de confiança para as médias.

alpha <- 0.05

star_estimates <- STAR %>%
  filter(!is.na(g4reading), classtype %in% classtypes_used) %>%
  group_by(classtype) %>%
  summarize(n = n(),
            est = mean(g4reading),
            se = sd(g4reading)/sqrt(n)) %>%
  mutate(lwr = est + qnorm(alpha/2)*se,
         upr = est + qnorm(1-alpha/2)*se)

star_estimates

## # A tibble: 2 × 6
##   classtype         n   est    se   lwr   upr
##   <fct>         <int> <dbl> <dbl> <dbl> <dbl>
## 1 Small class     726  723.  1.91  720.  727.
## 2 Regular class   836  720.  1.84  716.  723.

• Os intervalos possuem uma região em comum, isso sugere que é provável que o parâmetro dos processos geradores de dados seja o mesmo nas duas distribuições. A diferença observada na média pode ser por conta de características das amostras que estamos usando.
• Para investigar melhor a questão vamos calcular o intervalo de confiança para a diferença das médias.

star_ate <- star_estimates %>%
  arrange(desc(classtype)) %>% 
  summarize(se = sqrt(sum(se^2)),
            est = diff(est)) %>%
  mutate(ci_lwr = est + qnorm(alpha/2)*se,
         ci_up = est + qnorm(1-alpha/2)*se)

star_ate

## # A tibble: 1 × 4
##      se   est ci_lwr ci_up
##   <dbl> <dbl>  <dbl> <dbl>
## 1  2.65  3.50  -1.70  8.70

• O efeito médio do tratamento para classes pequenas foi estimado em 3,5 com erro padrão de 2,65. O intervalo de confiança de 95% é [-1,7, 8,7], que contém o zero.
• O resultado sugere que o efeito estimado, mesmo sendo positivo, tem grande chance de ser devido ao acaso. Não podemos afirmar que classes pequenas melhoram os resultados nos testes de leitura.

• Até agora calculamos intervalos de confiança usando o teorema central do limite. Como esse teorema se aplica a uma grande variedade de distribuições podemos construir intervalos de confiança usando distribuição normal sem muitas preocupações a respeito da verdadeira distribuição do processo gerador de dados.
• Ao usar o teorema central do limite estamos usando iferência assintótica. O problema com essa abordagem é que a amostra precisa ser suficientemente grande para fazer valer o teorema.

• Uma abordagem alternativa consiste em supor que a variável de resultado tem uma distribuição normal. Vamos usar o experimento STAR para ilustrar essa abordagem, para isso vamos supor que as notas dos testes de cada tipo de turma possuem distribuição normal.

• Uma olhada no histograma das notas não sugere que a distribuição normal seja uma boa hipótese. Ocorre que a inferência realizada com a abordagem de supor que a variável de resultado é normal é mais conservadora do que a abordagem assintótica, em nome desse conservadorismo alguns pesquisadores preferem supor que a distribuição é normal, mesmo quando não parece ser uma hipótese razoável, do que usar a abordagem assintótica.

• Quando a variável tem uma distribuição normal podemos usar a distribuição t de Student ou simplesmente \(t\) para construir intervalos de confiança.
• A notação \(t_v\) será usada para uma distribuição \(t\) com \(v\) graus de liberdade.

• Suponha que \(\{X_1, X_2, \dots, X_n\}\) são i.i.d. \(N(\mu, \sigma^2)\). Então o z-score da média da amostra, \(\bar{X}_n\), que é chamado de estatística-t, segue uma distribuição t de Student com \(n-1\) graus de liberdade. \[\mbox{estatística t da média da amostra}=\frac{\bar{X}_n-\mbox{média}}{\mbox{erro padrão}} \sim t_{n-1}\]

• A distribuição \(t\) é semelhante à de uma normal padrão, porém com caudas mais largas. A medida que os graus de liberdade aumentam, a distribuição \(t\) vai ficando mais parecida com a normal padrão.
• O valor esperado de uma distribuição \(t\) com \(v\) graus de liberdade é zero e a variância é \(\frac{v}{v-2}\) para \(v > 2\), se \(1 < v \leq 2\) a variãncia é \(\infty\) e, se \(v=1\), a variância é indefinida.

• A construção do intervalo de confiança segue os mesmos passos da abordagem assintótica, porém usando a distribuição \(t\) no lugar da normal. A mudança na distribuição altera os valores críticos usados nos intervalos.

• Valores critícos com 5 graus de liberdade

alfa <- c(0.01, 0.05, 0.1)
data.frame(alfa, normal = qnorm(1-alfa/2), t=qt(1-alfa/2, df=5))

##   alfa   normal        t
## 1 0.01 2.575829 4.032143
## 2 0.05 1.959964 2.570582
## 3 0.10 1.644854 2.015048

• Valores critícos com 50 graus de liberdade

data.frame(alfa, normal = qnorm(1-alfa/2), t=qt(1-alfa/2, df=50))

##   alfa   normal        t
## 1 0.01 2.575829 2.677793
## 2 0.05 1.959964 2.008559
## 3 0.10 1.644854 1.675905

• Construção dos intervalos de confiança usando a hipótese de que a variável de resultado tem distribuição normal.

alpha <- 0.05

star_estimates_t <- STAR %>%
  filter(!is.na(g4reading), classtype %in% classtypes_used) %>%
  group_by(classtype) %>%
  summarize(n=n(),
            est = mean(g4reading),
            se = sd(g4reading)/sqrt(n)) %>%
  mutate(lwr = est + qt(alpha/2, df=n-1)*se,
         upr = est + qt(1-alpha/2, df=n-1)*se)

star_estimates_t

## # A tibble: 2 × 6
##   classtype         n   est    se   lwr   upr
##   <fct>         <int> <dbl> <dbl> <dbl> <dbl>
## 1 Small class     726  723.  1.91  720.  727.
## 2 Regular class   836  720.  1.84  716.  723.

• O intervalo ficou muito próximo do obtido com o teorema do limite cenral (a diferença está nos decimais) porque temos muitos graus de liberdade.

• Para calcular o intervalo de confiança do estimador de diferenças de média suponha que \(\{X_1, X_2, \dots, X_n\}\) são i.i.d. \(N(\mu_X, \sigma^2_X)\) e \(\{Y_1, Y_2, \dots, Y_n\}\) são i.i.d. \(N(\mu_Y, \sigma^2_Y)\), a estatística \(t\) é dada por: \[\mbox{estatística t para diferença de médias}= \frac{(\bar{X}_n-\bar{Y}_m)-(\mu_X - \mu_Y)}{\sqrt{\frac{\sigma_X^2}{n}+\frac{\sigma_Y^2}{m}}}\]
• A estatística segue uma distribuição \(t\), mas o cálculo dos graus de liberdade é mais complicado e está além do escopo deste curso.

• De toda forma, podemos usar a função t.test() para construir o intervalo de confiança.

reading_small <- filter(STAR, classtype == "Small class")$g4reading
reading_reg <- filter(STAR, classtype == "Regular class")$g4reading

t.ci <- t.test(reading_small, reading_reg)

t.ci

## 
##  Welch Two Sample t-test
## 
## data:  reading_small and reading_reg
## t = 1.3195, df = 1541.2, p-value = 0.1872
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.703591  8.706055
## sample estimates:
## mean of x mean of y 
##  723.3912  719.8900

• Na Unidade 6 vimos o exemplo da nota de veto enviada por Schwarzenegger e queríamos saber se a mensagem formada pelas primeiras letras de cada linha podia ser fruto do acaso.
• Nossa estratégia foi calcular a probabilidade da combinação específica aparecer se as palavras que iniciam cada linha fossem escolhidas aleatoriamente dentre as palavras usadas na nota, ou seja, queríamos saber a probabilidade da mensagem decorrer do acaso.
• A estratégia de calcular a probabilidade que um evento observado ocorra a partir de um determinado modelo probabilístico é o princípio dos testes de hipóteses.

• A idéia de teste de hipótese foi apresentada por Ronald Fisher no livro The Design of Experiments, publicado em 1935, no livro o autor propõe o experimento da lady apreciadora de chá.
• Durante o chá na Universidade de Cambridge uma lady afirmou que o gosto do chá muda a depender se o chá era colocado no leite ou leite era colocado no chá. Para avaliar essa afirmação, Fischer propôs um experimento aleatório onde 8 xicáras idênticas foram preparadas, em 4, escolhidas de forma aleatória, o chá foi colocado antes do leite e, nas outras 4, o leite foi colocado antes do chá.

• A lady foi convidada para dizer como foi preparado o chá de cada xícara, especificamente ela tinha de separar as quatro xícaras com leite primeiro das quatro xícaras com chá primeiro. Ela acertou todas.
• Obra do acaso ou a lady realmente sabe diferenciar pelo sabor quando o chá é colocado antes do leite?

• Para analisar esse experimento vamos definir os resultados potenciais. Para cada uma das 8 xícaras, vamos considerar dois palpite potenciais da lady. Esses palpites podem ou não depender do leite ter sido colocado antes do chá.
• Com a hipótese que a lady não sabe diferenciar se primeiro foi o leite ou o chá, então os palpites dela não dependem do chá ou o leite terem sido colocados primeiro na xícara. Em outras palavras, com esta hipotese os dois resultado potenciais, chá -> leite e leite -> chá, são idênticos.

• O procedimento de Ficher foi calcular o número de palpites corretos em cada combinação possível com a hipótese que a lady não sabe reconhecer se primeiro veio o leite ou o chá.
• Na unidade anterior chamamos esse procedimento, onde o número de observações associados com cada condição é fixo, de aleatorização completa. Um experimento alternativo, onde a cada xícara seria sorteado quem colocar primeiro, leite ou chá, seria aleatorização simples.
• A próxima tabela ilustra o procedimento adotado por Fisher, cada cenário corresponde a uma possível ordem da xícaras com leite primeiro, M, ou chá primeiro, T.

• Como já tinha sido dito, a lady acertou todas as xícaras. No primeiro cenário teria acertado 4, no segundo 6, no terceiro 2 e assim por diante.
• Para calcular a probabilidade da lady ter acertado todas sob a hipótese que ela não sabe diferenciar os dois casos temos que calcular todas as combinações possíveis: \(_8C_4 = \frac{8!}{4! \times (8-4)!}=70\).
• A probabilidade da lady ter acertado por conta do acaso é de \(\frac{1}{70} \approx 0,\!01\), que é baixa. Parece que ela sabe mesmo reconhecer.

• Calculamos a probabilidade da lady acertar oito das oito xícaras, também podemos calcular a probabilidade dela acertar qualquer número entre 0 e 7. De saída, repare que a probabilidade de errar todas, acertar 0, é igual a probabilidade de acertar todas. Nos dois casos só existe uma combinação correta.
• Se entre as quatro apontadas com leite primeiro tiver uma certa e três erradas, no grupo das quatro apontadas com chá primeiro também tem uma certa e três erradas.

• De quantas maneiras isso pode acontecer? Pelo número de maneiras de acertar uma, \(_4C_1=\frac{4!}{1!(4-1)!}=4\), multiplicado pelo número de maneiras de errar três, \(_4C_3=\frac{4!}{3!(4-3)!}=4\), ou seja, o número de maneiras é 16.
• Pela mesma lógica o número de maneiras de acertar quatro (duas em cada grupo) é igual a \(_4C_2 \times _4C_2 = 36\). - Para acertar seis é preciso acertar três e errar uma de cada grupo,o número de maneiras que isso pode acontecer é 16.
• A probabilidade de cada número de acertos é igual ao número de maneiras que pode ocorrer, \(\{1, 16, 36, 16, 1\}\), dividido pelo número de combinações possíveis, 70.

cups <- 4
k <- c(0, seq_len(cups))

true <- tibble(correct = k*2,
               n = choose(cups, k)*choose(cups, cups-k)) %>%
  mutate(prob = n/sum(n))

true

## # A tibble: 5 × 3
##   correct     n   prob
##     <dbl> <dbl>  <dbl>
## 1       0     1 0.0143
## 2       2    16 0.229 
## 3       4    36 0.514 
## 4       6    16 0.229 
## 5       8     1 0.0143

• Pela natureza do experimento foi possível clacular as probabilidades analiticamente, nem sempre esse é o caso.
• Vamos usar Monte Carlo para fazer um cálculo aproximado das probabilidades no experimento do chá.
• Comecemos criando uma função para escolher as xícaras de forma aleatória e checar quantos palpites foram corretos.

sims <- 500
guess <- tibble(guess = c("M", "T", "T", "M", "M", "T", "T", "M"))

randomize_tea <- function(df) {
  assigment <- slice_sample(df,n=8) %>%
    rename(actual = guess)
  bind_cols(df, assigment) %>%
    summarize(correct = sum(guess==actual))
}

• A função slice_sample() sorteia \(n\) elementos do data.frame de forma aleatória, o deafult é sem reposição, caso queira com reposição é preciso usar o argumento replace como TRUE.

map_df(seq_len(10), ~ randomize_tea(guess))

## # A tibble: 10 × 1
##    correct
##      <int>
##  1       4
##  2       4
##  3       4
##  4       6
##  5       4
##  6       0
##  7       4
##  8       4
##  9       6
## 10       2

• Para nossa simulação vamos aplicar a função 1000 vezes e calcular a probabilidade pela frequência de cada número de acertos.

approx <- map_df(seq_len(sims), ~ randomize_tea(guess)) %>%
  count(correct) %>%
  mutate(prob = n/sum(n))

approx

## # A tibble: 5 × 3
##   correct     n  prob
##     <int> <int> <dbl>
## 1       0    11 0.022
## 2       2   107 0.214
## 3       4   254 0.508
## 4       6   122 0.244
## 5       8     6 0.012

• Agora vamos comparar o resultado aproximado, obtido com simulação, com o resultado exato, obtido analiticamente.

results <- approx %>%
  select(correct, prob_sim = prob) %>%
  left_join(select(true, correct, prob_exact = prob), by="correct") %>%
  mutate(diff = prob_sim - prob_exact)

results

## # A tibble: 5 × 4
##   correct prob_sim prob_exact     diff
##     <dbl>    <dbl>      <dbl>    <dbl>
## 1       0    0.022     0.0143  0.00771
## 2       2    0.214     0.229  -0.0146 
## 3       4    0.508     0.514  -0.00629
## 4       6    0.244     0.229   0.0154 
## 5       8    0.012     0.0143 -0.00229

• O exemplo da lady do chá nos ensina a lógica do teste de hipótese, uma espécie de versão probabilística da prova por contradição.
• Começamos o teste supondo com válida a hipótese que gostaríamos de refutar. Essa hipótese é chamada de hipótese nula e normalmente denotada por \(H_0\). No exemplo do chá, a hipótese nula é de que a lady não sabe diferenciar se primeiro colocaram leite o chá.

• Depois definimos uma estatística de teste, que é uma função dos dados observados. No exemplo, a estatísitica de teste é o número de classificações corretas. Na sequência, supondo a validade da hipótese nula, calculamos a distribuição da estatística de teste que por vezes é chamada de distribuição de referência.
• No exemplo calculamos a probabilidade de cada número de acerto supondo que a lady não sabe difenciar como o chá foi preparado.

• Finalmente, usamos a distribuição de referência para calcular a probabilidade do valor observado para a estatística. Assim podemos avaliar um quão provável é o valor observado caso a hipótese nula seja verdaderia. Se for provável não rejeitamos a hipótese nula, se for muito pouco provável rejeitamos a hipótese nula.

• Como quantificar “pouco provável”? Para isso usamos o p-valor.
• O p-valor pode ser entendido como a probabilidade de se observar um valor para estística de teste com módulo maior ou igual ao encontrado supondo que a hipótese nula é válida. Quanto menor o p-valor mais forte a evidência contrária à hipótese nula.
• Importante: o p-valor não representa a probabilidade da hipótese nula ser verdadeira.

• Para fazer um teste de hipótese deve especificar o nível do teste que será denotado, não por acaso, como \(\alpha\).
• Se o p-valor for menor ou igual a \(\alpha\) rejeitamos a hipótese nula, valores comuns para \(\alpha\) são 0,01 e 0,05.

• O nível do teste, \(\alpha\), representa a probabilidade de errar por rejeitar a hipótese nula quando esta hipótese é verdadeira, ou seja ocorre uma falsa rejeição de \(H_0\). Esse erro é chamado de erro do Tipo I. Quando não rejeitamos \(H_0\) mesmo com a hipótese nula sendo falsa temos o erro do Tipo II.

• Podemos especificar a probabilidade de cometer erro do Tipo I escolhendo o valor adequado para \(\alpha\), não temos como afetar diretamente o erro do Tipo II.
• Existe um trade-off entre os dois tipos de erro, via de regra minimizar o erro do Tipo I aumenta o erro do Tipo II.

• No exemplo do chá, a estatística de teste é o número de xícaras classificadas corretamente. Como o observado foi 8, o maior valor possível, o p-valor é igual a probabilidade de acertar 8 sob a hipótese nula, no caso \(\frac{1}{70} \approx 0,\!014\)
• Se a lady tivesse classificado corretamente 6 xícaras, teríamos dois valores pelo menos tão extremos (maiores) do que 6: 6 e 8. O p-valor seria: \(\frac{16}{70}+\frac{1}{70} \approx 0,\!243\).

• Esse tipo de p-valor que calculamos é chamado de unicaudal (um lado da dsitribuição) porque só considera os valores da estatísitica de teste que são maiores ou iguais aos valores observados.
• O p-valor bicaudal (dois lados da distribuição) considera valores extremos nas duas caudas da distribuição. Se a distribuição for simétrica, o p-valor bicaudal é o dobro do p-valor unicaudal.

• O tipo de teste proposto por Fisher no exemplo do chá é conhecido como teste exato de Fisher, o exato é porque o p-valor é calculado exatamente e não por aproximação. Outra caracterísitca importante desse tipo de teste é que a validade do teste é justificada com base no caráter aleatório na distribuição do tratamento.
• A função fisher.test() faz o teste exato de Fisher. O argumento da função é uma matriz de contigência \(2\times 2\) com linhas e colunas representando a distribuição do tratamento e a variável de resultado.

• A matriz em \(x\) é o caso com todas as classificações corretas e a matriz em \(y\) com 6 classficações corretas.

x <- matrix(c(4,0,0,4), byrow = TRUE, ncol=2, nrow=2)
y <- matrix(c(3,1,1,3), byrow = TRUE, ncol=2, nrow=2)
rownames(x) <- colnames(x) <- rownames(y) <- 
  colnames(y) <- c("M", "T")

x

##   M T
## M 4 0
## T 0 4

y

##   M T
## M 3 1
## T 1 3

• Teste unicaudal para \(x\)

fisher.test(x, alternative = "greater")

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value = 0.01429
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  2.003768      Inf
## sample estimates:
## odds ratio 
##        Inf

• Teste bicaudal para \(x\)

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value = 0.02857
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.339059      Inf
## sample estimates:
## odds ratio 
##        Inf

• Teste unicaudal para \(y\)

fisher.test(y)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  y
## p-value = 0.4857
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    0.2117329 621.9337505
## sample estimates:
## odds ratio 
##   6.408309

• Em geral, um teste de hipótese estatístico consiste nos seguintes passos:
  • Especificar uma hipótese nula e uma hipótese alternativa;
  • Escolher uma estaística de teste e o nível do teste, \(\alpha\);
  • Obter a distribuição de referência, ou seja, a distribuição da estatísitica do teste supondo a hipótese nula como verdadeira;
  • Calcular o p-valor, unicaudal ou bicaudal a depender da hipótese alternativa;
  • Rejeitar a hipótese nula se o p-valor for menor ou igual a \(\alpha\), caso contrário não rejeitar a hipótese nula.

• Os princípios acima podem ser usados para implementar vários tipos de testes. Vamos apresentar dois tipos de testes paraticularmente importantes em ciências sociais: teste de uma amostra e teste de duas amostras.
• Teste de uma amostra para médias são usados para avaliar a hipótese nula de que a média de uma população é igual a um determinado valor. Teste de duas amostras são usados para avaliar a hipótese nula de que as médias de duas populações são iguais.
• Esse últimos são usados para avaliar se a média de resultado no grupo tratado é igual a média no grupo de controle.

• Para ilustrar o teste de uma amostra retornemos ao exemplo das pesquisas eleitorais. Suponha que a hipótese nula é que exatamente metade dos eleitores apoiam Obama, ou seja, \(H_0: p= 0,\!5\) e a hipótese alternativa é que a taxa de apoio a Obama não é de metade dos eleitores, ou seja, \(H_1: p \neq 0,\!5\).
• Suponha que fizemos uma pesquisa com 1018 elitores, \(n=1018\), nessa amostra 550 afirmaram apoiar Obama e o resto não. Desta forma \(\bar{X}_n = \frac{550}{1018} = 0,\!54\).
• A proporção encontrada na pesquisa, \(0,\!54\) é diferente da que temos como hipótese, \(0,\!5\), a questão é saber se essa diferença é estatisticamente significativa.

• Pelo nosso procedimento de testes de hipótese o primeiro passo é definir as hipóteses nula e alternativa, isso já foi feito, \(H_0: p = 0,\!5\) e \(H_1: p \neq 0,\!5\).
• O próximo passo é escolher a estatística de teste e o nível de significância, a estatística é a média da amostra, \(\bar{X}_n\), e o nível de significânica é tal que \(\alpha = 0,\!5\).
• Para a distribuição de \(\bar{X}_n\) vamos usar o teorema central do limite, de forma que \(\bar{X}_n \sim N\left(0,\!5,\frac{0,\!5(1-0,\!5)}{1018}\right)\). Note que calculamos a média e a variância supondo que \(p = 0,\!5\), ou seja, que \(H_0\) é verdadeira.

• Como o teste é bicaudal (estamos testando \(p=0,\!5\) contra \(p < 0,\!5\) ou \(p > 0,\!5\)) devemos calcular o p-valor como a probabilidade que sob a hipotese nula observarmos um valor mais extremo do que o valor observado.
• Especificamente, queremos saber a probabilidade de sob a hipótese nula, \(\bar{X}_n > 0,\!54\) ou \(\bar{X}_n < 0,\!46\).

• A função pnorm() calcula a probabilidade que desejamos.

n <- 1018
x.bar <- 550/n
se <- sqrt(0.5 * 0.5/n)

upper <- pnorm(x.bar, mean = 0.5, sd = se, lower.tail = FALSE)
lower <- pnorm(0.5 - (x.bar - 0.5), mean = 0.5, sd=se)

upper+lower

## [1] 0.01016866

2*upper

## [1] 0.01016866

• No caso de teste unicaudal, \(H_0: p = 0.5\) e \(H_1: p >0.5\), o p-valor será:

upper

## [1] 0.005084332

• Nos dois testes, bicaudal e unicaudal, o p-valor foi menor do que 0,05, portanto, nos dois testes, rejeitamos a hipótese nula de \(p = 0,\!5\).
• O valor acima de 0,5 que encontramos não parece ter sido por conta do acaso.

• Quando usamos a distribuição normal como a distirbuição de referência, é comum usar o z-score para padronizar a estatística de teste, com essa transformação da distribuição se torna uma normal padrão. \[\frac{\bar{X}_n - \mu_0}{\mbox{erro padrão de }\bar{X}_n} \sim N(0,1)\]

• No exemplo:

z.score <- (x.bar-0.5)/se
z.score

## [1] 2.57004

pnorm(z.score, lower.tail = FALSE)

## [1] 0.005084332

2*pnorm(z.score, lower.tail = FALSE)

## [1] 0.01016866

• Suponha que \(\{X_1, X_2, \dots, X_N \}\) são variáveis aleatórias i.i.d com média \(\mu\) e variância \(\sigma^2\). O teste z de uma amostra consiste em:
  • Hipótese nula que a média populacional é igual a um valor especificado \(\mu_0\), \(H_0: \mu = \mu_0\).
  • Hipótese alternativa \(H_1: \mu \neq \mu_0\), teste bicaudal, \(H1: \mu > \mu_0\), teste unicaudal; ou \(H1: \mu < \mu_0\), teste unicaudal.
  • Estatística de teste (estatística z): \(Z_n = \frac{\bar{X}_n - \mu_0}{\sqrt{\frac{\hat{\sigma}^2}{n}}}\), onde \(\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i\).

• Continuação:
  • Para \(n\) grande a distribuição de referência é a \(N(0,1)\).
  • Variância: \(\hat{\sigma}^2=\frac{1}{n-1}\sum_{i=1}^n (X_i -\bar{X}_n)^2\) ou \(\bar{\sigma}^2=\frac{\mu_0(1-\mu_0}{n})\) se \(X_i\) for Bernoulli.
  • O p-valor será \(\Phi(-|Z_n|)\), unicaudal, ou \(2\Phi(-|Z_n|)\), bicaudal, onde \(\Phi()\) é a função de distribuição acumulada de uma normal padrão.

• Se \(X\) tem distribuição normal, a mesma estatística \(Z_n\) é chamada de estatística \(t\) e segue uma distribuição \(t\) com \(n-1\) graus de liberdade. O p-valor será calculado com a função de distribuição acumulada da distribuição \(t\). Esse teste é chamado teste t de uma amostra e é mais conservador do que o teste z.

• Existe uma relação próxima entre teste de hipótese e intervalo de confiança. Rejeitar a hipótese \(H_0: \mu = \mu_0\) em um teste bicaudal com nível de significância \(\alpha\) equivale a dizer que que \(\mu_0\) não pertence a u,m intervalo com nível de confiança \((1-\alpha)\times 100\%\).
• Para ilustrar esse resultado lembre que o p-valor do teste bicaudal do exemplo foi \(0,\!0102\), que é menor do que \(0,\!05\), porém maior do que \(0,\!01\), desta forma não rejeitamos \(H_0\) com 1% e rejeitamos com 5%, dito de outra forma, \(0,\!5\) deve pertencer ao intervalo com 99% de confiança e não pertencer ao intervalo com 95% de confiança.

c(x.bar - qnorm(0.995)*se, x.bar + qnorm(0.995)*se)

## [1] 0.4999093 0.5806408

c(x.bar - qnorm(0.975)*se, x.bar + qnorm(0.975)*se)

## [1] 0.5095605 0.5709896

• A função prop.test() implementa testes de hipótese e calcula o intervalo de confiança.

prop.test(550, n=n, p=0.5, correct = FALSE)

## 
##  1-sample proportions test without continuity correction
## 
## data:  550 out of n, null probability 0.5
## X-squared = 6.6051, df = 1, p-value = 0.01017
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5095661 0.5706812
## sample estimates:
##        p 
## 0.540275

• O argumento correct informa se deve ser feita uma correção para a proximação de uma distribuição binomial (discreta) por uma normal (contínua). O default do argumento é TRUE, em geral a correção é recomendada, especialmente quando \(n\) é pequeno.

prop.test(550, n=n, p=0.5)

## 
##  1-sample proportions test with continuity correction
## 
## data:  550 out of n, null probability 0.5
## X-squared = 6.445, df = 1, p-value = 0.01113
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5090744 0.5711680
## sample estimates:
##        p 
## 0.540275

• Podemos mudar o nível de significância (confiança) do teste (intervalo).

prop.test(550, n=n, p=0.5, conf.level = 0.99)

## 
##  1-sample proportions test with continuity correction
## 
## data:  550 out of n, null probability 0.5
## X-squared = 6.445, df = 1, p-value = 0.01113
## alternative hypothesis: true p is not equal to 0.5
## 99 percent confidence interval:
##  0.4994182 0.5806040
## sample estimates:
##        p 
## 0.540275

• No exemplo do projeto STAR tínhamos as notas dos estudantes no teste de leitura, suponha que queiramos testar se a média da população é igual a 710, \(H_0: \mu = 710\). Para realizar um teste \(t\) de uma amostra usamos a função t.test().

t.test(STAR$g4reading, mu=710)

## 
##  One Sample t-test
## 
## data:  STAR$g4reading
## t = 10.407, df = 2352, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 710
## 95 percent confidence interval:
##  719.1284 723.3671
## sample estimates:
## mean of x 
##  721.2478

• O outro tipo de teste que vamos discutir é o teste de duas amostras, esse teste é usado para, por exemplo, avaliar diferenças entre grupo tratado e grupo de controle.
• Vamos seguir com o exemplo do projeto STAR, a hipótese nula é que o efeito médio do tratamento na população, PATE, é zero, ou seja, \(H_0: E(Y_i(1)-Y_i(0))=0\), e a hipótese alternativa é que o PATE é diferente de zero, \(H_0: E(Y_i(1)-Y_i(0))\neq 0\).

• Caso existam fortes motivos para acreditar que o PATE não pode ser negativo, podemos usar um teste unicaudal com \(H_1: E(Y_i(1)-Y_i(0))>0\), da mesma forma, se temos fortes razões para acreditar que o PATE não pode ser positivo podemos fazer \(H_1: E(Y_i(1)-Y_i(0))<0\).
• No nosso exemplo vamos testar se o efeito de pequenas turmas é diferente de zero.

• A estatística de teste será o estimador de diferenças nas médias.
• Para distribuição de referência vamos usar o teorema central do limite, de forma que vamos considerar que as médias dos grupos de tratamento e controle tem distribuição normal.
• Pelas propriedades da normal a diferença das médias também segue uma distribuição normal. Com a hipótese nula, as médias são iguais nos dois grupos, a diferença das médias é normal com média zero. O z-score da diferença nas médias segue uma normal padrão.

• Suponha que \(\{X_1, X_2, \dots, X_{n_0}\}\) são variáveis aleatórias i.i.d com média \(\mu_0\) e variância \(\sigma^2_0\) e que \(\{Y_1, Y_2, \dots, Y_{n_1}\}\) são variáveis aleatórias i.i.d com média \(\mu_1\) e variância \(\sigma^2_1\). O teste z de duas amostras consiste em:
  • Hipótese nula que as duas populações têm a mesma média, \(H_0: \mu_0 - \mu_1\);
  • Hipótese alternativa \(H_1: \mu_0 \neq \mu_1\) bicaudal, \(H_1: \mu_0 > \mu_1\), unicaudal, ou \(H_1: \mu_0 < \mu_1\), unicaudal;
  • Estatística de teste: \(Z_n = \frac{\bar{Y}_{n_1}-\bar{X}_{n_0}}{\sqrt{\frac{1}{n_1}\hat{\sigma}_1^2+\frac{1}{n_0}\hat{\sigma}_0^2}}\);

• Continuação:
  • Distribuição de referência: \(Z_n \sim N(0,1)\) quando \(n_0\) e \(n_1\) são grandes o suficiente;
  • Variância: \(\hat{\sigma}_0^2 = \frac{1}{n_0-1}\sum_{i=1}^{n_0}(X_i - \bar{X}_{n_0})^2\) e \(\hat{\sigma}_1^2 = \frac{1}{n_1-1}\sum_{i=1}^{n_1}(Y_i - \bar{Y}_{n_1})^2\) ou \(\hat{\sigma}_0^2 = \hat{\sigma}_1^2=\hat{p}(1-\hat{p})\) com \(\hat{p}=\frac{n_0}{n_0+n_1}\bar{X}_{n_0}+\frac{n_1}{n_0+n_1}\bar{Y}_{n_1}\) se \(X\) e \(Y\) são Bernoulli;
  • O p-valor será \(\Phi(-|Z_n|)\), unicaudal, ou \(2\Phi(-|Z_n|)\), bicaudal, onde \(\Phi()\) é a função de distribuição acumulada de uma normal padrão.

• Se \(X\) e \(Y\) têm distribuição normal, a estatística \(Z_n\) é chamada estatística \(t\) e segue uma \(t\) de Student, nesse caso o p-valor será dado pela função de distribuição acumulada da distribuição \(t\).
• Esse teste é chamado teste t de duas amostras e é mais conservador do que o teste z de duas amostras.

• O data.frame star_ate tem o erro padrão, o valor estimado e os limites superior e inferior do intervalo de confiança.

head(star_ate)

## # A tibble: 1 × 4
##      se   est ci_lwr ci_up
##   <dbl> <dbl>  <dbl> <dbl>
## 1  2.65  3.50  -1.70  8.70

• Vamos acrescentar o p-valor para o teste z de duas amostras

star_ate %>%
  mutate(p_value_1sided = pnorm(-abs(est), mean=0, sd=se),
         p_value_2sided = 2*pnorm(-abs(est), mean=0, sd=se))

## # A tibble: 1 × 6
##      se   est ci_lwr ci_up p_value_1sided p_value_2sided
##   <dbl> <dbl>  <dbl> <dbl>          <dbl>          <dbl>
## 1  2.65  3.50  -1.70  8.70         0.0935          0.187

• Nos dois testes o p-valor ficou bem acima de 5%, desta forma, não rejeitamos a hipótese nula de que as médias dos dois grupos são iguais.
• O resultado não é uma surpresa, pois vimos que o intervalo de confiança contém o zero.

• A função t.test() pode ser usada para fazer o teste de t de duas amostras.

reading_small <- filter(STAR, classtype == "Small class")$g4reading
reading_reg <- filter(STAR, classtype == "Regular class")$g4reading

t.test(reading_small, reading_reg)

## 
##  Welch Two Sample t-test
## 
## data:  reading_small and reading_reg
## t = 1.3195, df = 1541.2, p-value = 0.1872
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.703591  8.706055
## sample estimates:
## mean of x mean of y 
##  723.3912  719.8900

• Para outro exemplo de teste de hipótese, retornemos ao caso da discriminação no mercado de trabalho discutido na Unidade 2.
• Queremos testar a hipótese de que a probabilidade de receber uma ligação de retorno é a mesma para nomes associados a afro-americanos e nomes associados a americanos brancos.

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

##   firstname    sex  race call
## 1   Allison female white    0
## 2   Kristen female white    0
## 3   Lakisha female black    0
## 4   Latonya female black    0
## 5    Carrie female white    0
## 6       Jay   male white    0

x <- resume %>%
  count(race, call) %>%
  pivot_wider(names_from = call, values_from = n) %>%
  ungroup()

x

## # A tibble: 2 × 3
##   race    `0`   `1`
##   <chr> <int> <int>
## 1 black  2278   157
## 2 white  2200   235

prop.test(as.matrix(select(x,-race)), alternative = "greater")

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  as.matrix(select(x, -race))
## X-squared = 16.449, df = 1, p-value = 2.499e-05
## alternative hypothesis: greater
## 95 percent confidence interval:
##  0.01881967 1.00000000
## sample estimates:
##    prop 1    prop 2 
## 0.9355236 0.9034908

• O p-valor é bem menor que 0,05, também é menor que 0,01, ou seja, rejeitamos a hipótese nula que a probabilidade de receber um retorno é a mesma para nomes associados a brancos e afro-americanos.

• Já fizemos o teste, mas, ainda assim, é válido calcular o p-valor passo a passo. Com a hipótese nula temos \(\mu_0=\mu_1\) e o erro padrão é dado por: \[\sqrt{\frac{\widehat{V(X)}}{n_0}+\frac{\widehat{V(Y)}}{n_1}}=\sqrt{\frac{\hat{p}(1-\hat{p})}{n_0}+\frac{\hat{p}(1-\hat{p})}{n_1}}=\\\sqrt{\hat{p}(1-\hat{p})\left(\frac{1}{n_0}+\frac{1}{n_1}\right)}\] \[\hat{p}=\frac{\sum_{i=1}^{n_0}X_i + \sum_{i=1}^{n_1}Y_i}{n_0 + n_1}\]

n0 <- sum(resume$race == "black")
n1 <- sum(resume$race == "white")

p <- mean(resume$call)
p0 <- mean(filter(resume, race == "black")$call)
p1 <- mean(filter(resume, race == "white")$call)

est <- p1 - p0
est

## [1] 0.03203285

se <- sqrt(p*(1-p)*(1/n0 + 1/n1))
se

## [1] 0.007796894

zstat <- est/se
zstat

## [1] 4.108412

pnorm(-abs(zstat))

## [1] 1.991943e-05

• O p-valor não ficou idêntico por conta da correção realizada pela função prop.test(), podemos resolver isso.

prop.test(as.matrix(select(x,-race)), alternative = "greater", correct = FALSE)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  as.matrix(select(x, -race))
## X-squared = 16.879, df = 1, p-value = 1.992e-05
## alternative hypothesis: greater
## 95 percent confidence interval:
##  0.01923035 1.00000000
## sample estimates:
##    prop 1    prop 2 
## 0.9355236 0.9034908

• O problema dos testes múltiplos está relacioanado a questão de que quando conduzimos muitos testes de hipóteses é provável encontrarmos conclusões falsas, especificamente, é provável rejeitarmos hipóteses nulas que são verdadeiras.

• Na Copa de 2010 o polvo Paul fez sucesso prevendo resultados dos jogos da Alemanha.
• Paul acertou a o vencedor dos sete jogos disputados pela Alemanha e acertou o vencedor da final da Copa.

• Podemos testar os poderes de Paul. Para isso a hipótese nula é de que Paul acertou na sorte, qual a probabilidade de acertar os oito resultados sob a hipótese nula que Paul não tem poderes especiais?
• Se o resultado for apenas vitória ou derrota temos uma binomial com probabilidade de sucesso 0,5 e tamanho 8.

dbinom(8, size=8, prob = 0.5)

## [1] 0.00390625

1/(2^8)

## [1] 0.00390625

• Com um pouco mais de preciosismo podemos lembrar que nos três jogos da primeira fase é possível ocorrer empate, nesse caso a probabilidade é dada por:

dbinom(3, size=3, prob = 1/3) * dbinom(5, size=5, prob=0.5)

## [1] 0.001157407

1/(3^3) * 1/(2^5)

## [1] 0.001157407

• Nos dois casos o p-valor ficou abaixo de 0,05, o que nos levaria a rejeitar a hipótese nula de que Paul não tem poderes especiais. Será?
• Qual a chance de rejeitar \(H_0\) mesmo ela estando correta, ou seja, qual a chance de cometer um erro do Tipo I? A resposta é \(\alpha\), que em geral é 5%. Como a probabilidade de rejeitar \(H_0\) mesmo sendo verdadeira é \(\alpha\), a probabilidade de não rejeitar \(H_0\) quando ela é verdadeira é \(1-\alpha\), que em geral é 95%.

• A probabilidade de não rejeitar \(H_0\) quando ela é verdadeira duas vezes em dois testes é \(0,\!95^2\). Seguindo com essa lógica a probabilidade de não rejeitar \(H_0\) quando ela é verdadeira dez vezes em dez testes é de \(0,\!95^{10} \approx 0,\!6\).
• Portanto, a probabilidade de rejeitar \(H_0\) pelo menos uma vez em dez testes mesmo ela sendo verdadeira é de \(1-0,\!95^{10} \approx 0,\!4\). Talvez não seja tão difícil encontrar outro Paul…

• O poder de um teste de hipótese é definido como a probabilidade de rejeitar a hipótese nula quando a hipótese nula é falsa, ou seja, o poder de um teste é igual a um menos a probabilidade do erro de Tipo II. \[\mbox{poder}=1-P(\mbox{erro do Tipo II})\]

• O poder do teste permite analisar o grau de informação dos dados usados para fazer o teste de hipótese.
• Em particular é usado para definir o menor tamanho da amostra necessária para estimar o parâmetro com precisão suficente para que seja possível distinguir o valor observado da estatística do valor suposto na hipótese nula.

• Voltemos ao exemplo da pesquisa eleitoral. Suponha que queiramos determinar quantas pessoas devem ser entrevistadas para sermos capazes de rejeitar a hipótese nula que Obama tem exatamente 50% de apoio, \(p=0,\!5\), quando o apoio verdadeiro está a pelo menos dois pontos percentuais de \(0,\!5\), ou seja, está entre \(0,\!48\) e \(0,\!52\).
• Dito de outra forma, dois pontos percentuais é o menor desvio da hipótese nula que queremos detectar com alta probabilidade.
• Considere que a estatística de teste é a proporção encontrada na amostra, que o nível de significância é \(\alpha = 0,\!05\) e o teste é bicaudal.

• Para calcular o poder do teste devemos considerar a distribuição sob a hipótese nula e a distribuição com a proporção observada na amostra. Seja \(p\) a proporção sob a hipótese nula e \(p^\star\) a proporção obervada na amostra, então as distribuições são: \(N\left(p, \frac{p(1-p)}{n} \right)\) e \(N\left(p^\star, \frac{p^\star(1-p^\star)}{n} \right)\).
• Suponha \(p=0,\!5\), \(p^\star=0,\!48\) e \(n=250\).

• A distribuição em preto foi construída com \(p\), a distribuição em azul com \(p^\star\). As linhas pontilhadas representam os valores extremos para rejeição de \(H_0\).
• O poder do teste é a soma das áreas destacadas, a esquerda da fronteira inferior e a direita da superior abaixo da linha azul. Formalmente: \[\mbox{poder}= P(\bar{X}_n \leq p-z_{\alpha/2}\times\mbox{erro padrão}) \\+ P(\bar{X}_n \geq p+z_{\alpha/2}\times\mbox{erro padrão})\]

• Vamos usar o R para calcular o poder do teste

n <- 250
p.star <- 0.48
p <- 0.5
alpha <- 0.05

cr.value <- qnorm(1-alpha/2)
se.star <- sqrt(p.star * (1-p.star)/n)
se <- sqrt(p*(1-p)/n)

#poder do teste
pnorm(p - cr.value * se, mean = p.star, sd=se.star) + 
  pnorm(p + cr.value * se, mean = p.star, sd=se.star, lower.tail = FALSE)

## [1] 0.09673114

• O poder do teste depende do tamanho da amostra:

n <- 5000
p.star <- 0.48
p <- 0.5
alpha <- 0.05

cr.value <- qnorm(1-alpha/2)
se.star <- sqrt(p.star * (1-p.star)/n)
se <- sqrt(p*(1-p)/n)

#poder do teste
pnorm(p - cr.value * se, mean = p.star, sd=se.star) + 
  pnorm(p + cr.value * se, mean = p.star, sd=se.star, lower.tail = FALSE)

## [1] 0.8076207

• O poder do teste depende de \(p\):

n <- 250
p.star <- 0.48
p <- 0.3
alpha <- 0.05

cr.value <- qnorm(1-alpha/2)
se.star <- sqrt(p.star * (1-p.star)/n)
se <- sqrt(p*(1-p)/n)

#poder do teste
pnorm(p - cr.value * se, mean = p.star, sd=se.star) + 
  pnorm(p + cr.value * se, mean = p.star, sd=se.star, lower.tail = FALSE)

## [1] 0.9999517

• O poder do teste é definido como a probabilidade de rejeitar a hipótese nula quando a hipótese nula é falsa, que é igual a um menos a probabilidade do erro do Tipo II. Análise de poder do teste consiste em:
  • Selecionar as hipóteses estatísticas que serão usadas: especificação da estatística de teste, hipóteses nula e alternativa, nível de signficânica, etc;
  • Escolher o parâmetro populacional sob um processo gerador de dados hipotético;
  • Calcular a probabilidade de rejeitar a hipótese nula sob o suposto processo gerador de dados e um determinado tamanho da amostra.
• Então escolhemos o tamanho da amostra necessário para o pode do teste desejado.

• A análise de poder pode ser feita em testes de duas amostras. Considere que você aplicou um tratamento em um grupo e deseja comparar proporções no grupo tratado com a de um grupo de controle.
• A estatística de teste é a diferença entre as proporções dos grupos de tratamento e controle, \(\bar{Y}_{n_1}-\bar{X}_{n_0}\). Sob a hipótese nula a diferença entre as proporções da duas populações, PATE, é zero.

• A distribuição de referência é \(N\left(0, p(1-p)\left(\frac{1}{n_1}+\frac{1}{n_0}\right)\right)\) onde \(p=\frac{n_0 p_0 + n_1p_1}{n_0 + n_1}\).
• Para calcular o poder do teste devemos especificar as proproções separadamente para cada grupo, desta forma a distribuição será: \(N\left(p_1^\star-p_0^\star, \frac{p_1^\star(1-p_1^\star)}{n_1}+\frac{p_0^\star(1-p_0^\star)}{n_0}\right)\).

• Como exemplo considere o caso dos currículos. Suponha que planejamos mandar 500 currículos com nomes associados a afro-americanos e outros 500 com nomes associados a brancos. Suponha também que esperamos 5% de retorno para os afro-americanos e 10% para os brancos.

n1 <- 500
n0 <- 500
p1.star <- 0.05
p0.star <- 0.1

• Para calcular o poder do teste primeiro calculamos o retorno médio, usando a média ponderada dos dois grupos, \(p = \frac{n_1p_1^\star + n_0p_0^\star}{n_1+n_0}\).
• Depois calculamos o erro padrão sob a hipótese nula, \(\mbox{se}=\sqrt{p(1-p)\left(\frac{1}{n_1}+\frac{1}{n_0}\right)}\).
• Finalmente calculamos o erro padrão com \(p_1^\star\) e \(p_2^\star\), \(\mbox{se}^\star=\sqrt{\frac{p_1^\star(1-p_1^\star)}{n_1}+\frac{p_0^\star(1-p_0^\star)}{n_0}}\).

p <- (n1*p1.star + n0*p0.star)/(n1+n0)
se <- sqrt(p*(1-p)*(1/n1 + 1/n0))
se.star <- sqrt(p1.star*(1-p1.star)/n1 + p0.star*(1-p0.star)/n0)

• Poder do teste:

pnorm(-cr.value*se, mean = p1.star - p0.star, sd=se.star) +
  pnorm(cr.value*se, mean=p1.star-p0.star, sd=se.star, lower.tail = FALSE)

## [1] 0.85228

• A função power.prop.test() calcula o poder do teste de proporções com duas amostras.
• A função também calcula a tamanho da amostra para um dado poder do teste.

power.prop.test(n=500, p1=0.05, p2=0.1, sig.level = 0.05)

## 
##      Two-sample comparison of proportions power calculation 
## 
##               n = 500
##              p1 = 0.05
##              p2 = 0.1
##       sig.level = 0.05
##           power = 0.8522797
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

power.prop.test(p1=0.05, p2=0.1, sig.level = 0.05, power=0.9)

## 
##      Two-sample comparison of proportions power calculation 
## 
##               n = 581.0821
##              p1 = 0.05
##              p2 = 0.1
##       sig.level = 0.05
##           power = 0.9
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

• Para variáveis aleatórias contínuas a função power.t.tes() calcula poder de teste tanto de uma como de duas amostras.
• No teste de uma amostra o argumento type deve ser definido como one.sample, o argumento delta como a média e o argumento sd como o desvio padrão de um hipótetico processo gerador de dados.
• Para o teste de duas amostras é preciso que o tamanho da amostra e o desvio padrão são iguais nos dois grupos. O argumento delta especifica a diferença nas médias.

• Calcular o poder de um teste de uma amostra de tamanho 100, com média igual a 0,25 e desvio padrão igual a 1, na hipótese nula a média é igual a zero.

power.t.test(n=100, delta = 0.25, sd=1, type="one.sample")

## 
##      One-sample t test power calculation 
## 
##               n = 100
##           delta = 0.25
##              sd = 1
##       sig.level = 0.05
##           power = 0.6969757
##     alternative = two.sided

• No teste acima, qual o tamanho da amostra para o poder ser 90%?

power.t.test(power=0.9, delta = 0.25, sd=1, type="one.sample")

## 
##      One-sample t test power calculation 
## 
##               n = 170.0511
##           delta = 0.25
##              sd = 1
##       sig.level = 0.05
##           power = 0.9
##     alternative = two.sided

• Duas amostras com diferença de média igual a 0,25.

power.t.test(n=100, delta = 0.25, sd=1, type="two.sample")

## 
##      Two-sample t test power calculation 
## 
##               n = 100
##           delta = 0.25
##              sd = 1
##       sig.level = 0.05
##           power = 0.4204383
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

• No teste acima, qual o tamanho da amostra para o poder ser 90%?

power.t.test(power=0.9, delta = 0.25, sd=1, type="two.sample")

## 
##      Two-sample t test power calculation 
## 
##               n = 337.2008
##           delta = 0.25
##              sd = 1
##       sig.level = 0.05
##           power = 0.9
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

• Para terminar o curso voltemos ao modelo de regressão linear. Porém, desta vez vamos analisar o modelo como uma aproximação para o processo gerador de dados. Com essa abordagem vamos quantificar a incerteza de nossas estimativas.

• O modelo de regressão linear com \(p\) variáveis explicativas (que podem ser chamadas de independentes ou preditoras) é definido por: \[Y_i = \alpha + \beta_1 X_{i1}+\beta_2X_{i2}+\cdots+\beta_pX_{ip} + \varepsilon_i\]
• No modelo temos \(p+1\) coeficientes que devemos estimar.

• De acordo com o modelo a variável de resultado (ou dependente) é gerada como uma função linear das variáveis explicativas e do termo de erro. Tanto a variável de resultado quanto as variáveis explicativas são observadas, apenas o erro não é observado.
• Por isso é fundamental fazer hipóteses a respeito da distribuição do erro.

• A hipótese de exogeneidade para o modelo de regressão linear é definida como: \[E(\varepsilon_i|X_1, X_2, \dots, X_n)=E(\varepsilon_i)=0\]
• Esta hipótese implica que a parte não observada do resultado, contida no termo de erro, não é relacionada com as variáveis explicativas \(X_{ij}\) para \(i=1,2,\dots,n\) e \(j=1,2,\dots,p\).

• De acordo com a hipótese de exogeneidade a epserança do erro condicional às variáveis explicativas é igual a esperança não condicional do erro.
• Desde que o modelo tenha uma constante, a esperança não condicional do erro pode sempre ser suposta igual a zero sem maiores problemas.

• A esperança condicional de uma variável aleatória \(Y\) dada outra variável aleatória \(X\) é denotada por \(E(Y|X)\) e dada por: \[E(Y|X)=\begin{cases}\sum_y y\times f(y|X) \qquad \ \mbox{se }Y\mbox{é discreta}\\ \int y\times f(y|X)dy\qquad \mbox{se }Y\mbox{é contínua}\end{cases}\] Onde \(f(y|X)\) é função de probabilidade condicional (função de densidade condicional) de uma variável aleatória \(Y\) discreta (contínua) dado \(x\).

• Com a hipótese de exogenidade temos que: \[E(Y_i|X_1,X_2,\dots,X_p)= \alpha + \beta_1 X_{i1}+\beta_2 X_{i2} + \cdots + \beta_p X_{ip}\]
• Para observar esse fato lembre que \(E(\beta_j X_{ij}|X_1, X_2, \dots, X_p)=\beta_j X_{ij}\)

• Em estudos aleatórios controlados, a hipótese de exogeneidade é respeitada porque a escolha das unidades tratadas é feita de forma aleatória. Em termos de modelo de regressão linear isso significa que a variável de tratamento, \(X\), é independente de todas as características, observadas ou não, que existem antes do tratamento que estão contidas em \(\varepsilon\).

• Considere o caso de mulheres e políticas públicas que vimos na Unidade 4, a variável explicativa de interesse, \(X\), era se existia reserva para mulheres no Gram Panchayat (GP) da vila. Como a distribuição das reservas foi aleatória, a vila ter a reserva não está relacionado a outras variáveis que podiam influenciar a variável de resultado como tamanho da população ou nível de renda.

• Em estudos com dados observados a validade da hipótese de exogeneidade é bem mais complicada. No exemplo dos GPs se não existisse o programa de reserva de vagas e usássemos a presença de mulheres no GP a exogeneidade poderia ficar comprometida.
• Imagine, por exemplo, que vilas com maior nível de educação tenham uma visão mais liberal e elejam mais mulheres e dêem mais valor a água tratada. Como saber se o maior investimento em equipamentos para taratmento da água foi por conta da presença de mulheres no GP ou pelas características dos eleitores?

• No estudos com dados observados, confundidores não observados (e.g. nível de educação) podem estar contidos no termo de erro, se também estiverem correlacionados com as variáveis explicativas (e.g. gênero dos políticos eleitos) a hipótese de exogeneidade fica comprometida.

• Existem várias técnicas para contornar o problema da não exogeneidade em estudos com dados observados. Uma delas, que discutimos na Unidade 2, é encontar unidades muito parecidas com as tratadas e que não receberam o tratamento.
• Foi a estratégia usada no estudo do salário mínimo em NJ, os autores do estudo usaram restaurantes das mesmas cadeias localizados na vizinha PA como controle para avaliar os efeito do aumento do salário mínimo em NJ.

• Em um modelo de regressão linear, se todos os confundidores estiverem no modelo (e a hipótese de relação linear se aplicar) o coeficiente estimado para a variável de tratamento é um estimador não viesado do efeito médio do tratamento.
• No exemplo do salário mínimo, suponha que os únicos confundidores entre os restaurantes fast-food em NJ e PA sejam a cadeia de fast-food, o salário e a proporção de empregados de tempo integral antes da mudança do salário mínimo em NJ.
• Vamos estimar o modelo com o emprego tempo integral como variável de resultado, estar em NJ como tratamento e as três variáveis acima como confundidores.

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

minwage <- minwage %>%
  mutate(fullPropBefore = fullBefore/(fullBefore + partBefore),
         fullPropAfter = fullAfter/(fullAfter + partAfter),
         NJ = if_else(location == "PA", 0, 1))

• Para que o R estime os coeficientes de cada cadeia de fast-food vamos retirar a constante do modelo.

fit_minwage <- lm(fullPropAfter ~ -1 + NJ + fullPropBefore + wageBefore + chain, 
                  data = minwage)

fit_minwage

## 
## Call:
## lm(formula = fullPropAfter ~ -1 + NJ + fullPropBefore + wageBefore + 
##     chain, data = minwage)
## 
## Coefficients:
##              NJ   fullPropBefore       wageBefore  chainburgerking  
##         0.05422          0.16879          0.08133         -0.11563  
##        chainkfc        chainroys      chainwendys  
##        -0.15080         -0.20639         -0.22013

• Se deixarmos a constante, o coeficiente dacadeia de referência, Burger King, é retirado e os coeficientes de cada cadeia passam a representar a diferença para o coeficiente do Burger King.

fit_minwage1 <- lm(fullPropAfter ~ NJ + fullPropBefore + wageBefore + chain, 
                   data = minwage)

fit_minwage1

## 
## Call:
## lm(formula = fullPropAfter ~ NJ + fullPropBefore + wageBefore + 
##     chain, data = minwage)
## 
## Coefficients:
##    (Intercept)              NJ  fullPropBefore      wageBefore        chainkfc  
##       -0.11563         0.05422         0.16879         0.08133        -0.03517  
##      chainroys     chainwendys  
##       -0.09076        -0.10451

• Os modelos estimados com e sem a constante são os mesmos, só muda a forma de apresentar os resultados. Repare, por exemplo, que o coeficiente da variável de interesse, \(NJ\), é o mesmo.
• Para ilustrar que os modelos são os mesmos vamos comparar as previsões de cada modelo.

library(modelr)

pred_1 <- minwage %>%
  slice(1) %>% 
  add_predictions(fit_minwage) %>%
  select(pred) %>%
  mutate(model = "fit_minwage")

pred_compare <- minwage %>%
  slice(1) %>%
  add_predictions(fit_minwage1) %>%
  select(pred) %>%
  mutate(model = "fit_minwage1") %>%
  bind_rows(pred_1)

pred_compare

##        pred        model
## 1 0.2709367 fit_minwage1
## 2 0.2709367  fit_minwage

• O uso do modelo linear para fazer inferência necessita da hipótese de exogeneidade. Essa hipótese será violada quando existem confundidores não observados. Para tornar a hipótese de exogeneidade mais plausível, pesquisadores medem variáveis confundidoras e as incluem como variáveis explicativas do modelo de regressão linear.

• Quão precisos são os coeficientes estimados no modelo de regressão linear? Com a hipótese que o modelo linear descreve o verdadeiro processo gerador de dados, podemos quantificar a incerteza associada aos coeficientes estimados.
• Por simplicidade, a análise será feita com a hipótese de apenas uma variável explicativa, mas os resultados valem para modelos com mais de uma variável explicativa.

• Modelo: \[Y_i = \alpha + \beta X_i + \varepsilon_i\]

• Coeficientes estimados \[\begin{align*}\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}\end{align*}\]

• Com a hipótese de exogeneidade \(\hat{\alpha}\) e \(\hat{\beta}\) são estimadores não viesados de \(\alpha\) e \(\beta\), ou seja, \(E(\hat{\alpha})=\alpha\) e \(E(\hat{\beta})=\beta\).

• Para mostrar esse resultado considere que: \[\begin{align*}\bar{Y} &= \alpha + \beta \bar{X} + \bar{\varepsilon} \Rightarrow \\ \hat{\alpha} &= \alpha + \beta \bar{X}+\bar{\varepsilon}-\hat{\beta}\bar{X}\Rightarrow \\ \hat{\alpha}&=\alpha +(\beta - \hat{\beta})\bar{X}+\bar{\varepsilon}\end{align*}\]
• O erro de estimação, \(\hat{\alpha}-\alpha\) é dado por \((\hat{\beta}-\beta)\bar{X}+\bar{\varepsilon}\)

• Para encontrar \(\hat{\beta}-\beta\):

\[\begin{align*}\hat{\beta}&=\frac{\sum_{i=1}^n (\beta X_i+\varepsilon_i-\beta\bar{X}-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow\\ \hat{\beta}&=\frac{\sum_{i=1}^n[\beta(X_i-\bar{X})+(\varepsilon_i-\bar{\varepsilon})](X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow\\ \hat{\beta}&=\frac{\sum_{i=1}^n \beta (X_i-\bar{X})^2}{\sum_{i=1}^n(X_i - \bar{X})^2}+\frac{\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow \\ \hat{\beta} &= \beta + \underbrace{\frac{\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}}_{\mbox{erro de estimação}}\end{align*}\]

• O numerador do erro de estimação acima pode ser escrito como:

\[\begin{align*}\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X}) &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X})+\sum_{i=1}^n\bar{\varepsilon} (X_i - \bar{X})\\ &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X}) + \bar{\varepsilon}\left(\sum_{i=1}^nX_i - n\bar{X} \right) \\ &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\end{align*}\]

• Desta forma: \[\hat{\beta}-\beta = \frac{\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})}{\sum_{i=1}^n(X_1-\bar{X})^2}\]

• Para seguir com o argumento vamos calcular a esperança condicional de \(\hat{\beta}-\beta\): \[E(\hat{\beta}-\beta|X)=\frac{1}{\sum_{i=1}^n(X_i - \bar{X})^2}\sum_{i=1}^nE(\varepsilon_i|X)(X_i-\bar{X})\]
• Pela hipótese de exogeneidade, \(E(\varepsilon_i|X) = 0\), logo: \[\begin{align*}E(\hat{\beta}-\beta|X) &= 0\\E(\hat{\alpha}-\alpha|X) &= E(\hat{\beta}-\beta|X)\bar{X}+E(\bar{\varepsilon}|X)=0\end{align*}\]
• Pois: \[E(\bar{\varepsilon}|X)=\frac{1}{n}\sum_{i=1}^nE(\varepsilon_i|X)=0\]

• A lei das expectativas iteradas diz que para quaisquer duas variáveis aleatórias \(X\) e \(Y\) vale que: \[E(Y)=E(E(Y|X))\]

• Pela lei das expectativas iteradas:

\[\begin{align*} E(\hat{\alpha})&=E(E(\hat{\alpha}|X))=E(\alpha)=\alpha\\E(\hat{\beta})&=E(E(\hat{\beta}|X))=E(\beta)=\beta\end{align*}\]

• Já temos os valores esperados para os estimadores, precisamos agora do erro padrão.

\[V(\hat{\beta}|X)=V(\hat{\beta}-\beta|X)=V\left(\frac{\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})}{\sum_{i=1}^n(X_i-\bar{X})^2}\bigg|X\right)\\=\frac{1}{\left[\sum_{i=1}^n(X_i - \bar{X})^2\right]^2}V\left(\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\bigg|X\right)\]

• Para simplificar a expressão acima usamos a hipótese de homocedasticidade, qual seja, condicional em \(X\) os erros \(\varepsilon_i\) são independentes e que a variância do erro não depende de \(X\).

• A hipótese de homocedasticidade dos erros consiste em:
  • \(\varepsilon_i\) é independente de \(\varepsilon_j\) condicional em \(X\) para todo \(i \neq j\);
  • A variância do erro não depende de \(X\), \(V(\varepsilon_i|X)=(\varepsilon_i)\)

• Com a hipótese de homocedasticidade, temos que:

\[V\left(\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\bigg|X\right)=\sum_{i=1}^n V(\varepsilon_i|X)(X_i-\bar{X})^2\\=V(\varepsilon_i)\sum_{i=1}^n(X_i-\bar{X})^2\] - Portanto: \[V(\hat{\beta}|X)=\frac{V(\varepsilon_i)}{\sum_{i=1}^n(X_i - \bar{X})^2}\]

• A lei da variância total diz que para quaisquer variáveis aleatórias \(X\) e \(Y\) vale que: \[V(Y)=V[E(Y|X)]+E[V(Y|X)]\]

• Aplicando a lei da variância total: \[V(\hat{\beta})=V[E(\hat{\beta}|X)]+E[V(\hat{\beta}|X)]=V(\beta)+E\left(\frac{V(\varepsilon_i)}{\sum_{i=1}^n(X_i - \bar{X})^2}\right)\\=V(\varepsilon_i)E\left(\frac{1}{\sum_{i=1}^n(X_i - \bar{X})^2}\right)\]

• Repare que o variância não condicional de \(\hat{\beta}\) é igual ao valor esperado da variância condicional de \(\hat{\beta}\), por conta disso um bom estimador da variância condicional também é um bom estimador da variância não condicional.

• Com esse resultado podemos calcular o erro padrão de \(\hat{\beta}\) com a raiz do estimador de variância. \[\mbox{erro padrão de }\hat{\beta}=\sqrt{\widehat{V(\hat{\beta})}}=\sqrt{\frac{\frac{1}{n}\sum_{i=1}^n \hat{\varepsilon}_i^2}{\sum_{i=1}^n(X_i - \bar{X})^2}}\]
• Para obter esse resultado usamos que média amostral dos resíduos é zero, de forma que: \(\frac{1}{n}\sum_{i=1}^n(\hat{\varepsilon}_i-\bar{\hat{\varepsilon}})^2=\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}_i^2\).

• Sem a hipótese de homocedasticidade a fórmula acima não é correta para estimar o erro padrão do estimador.
• Não é objetivo desse curso, mas o R tem pacotes com funções que calculam o erro padrão robusto para heterocedasticidade.

• Com o estimador para o erro padrão podemos construir intervalos de confiança para \(\hat{\beta}\). Usando o teorema do limite central temos que: \[\mbox{z-score de }\hat{\beta}=\frac{\hat{\beta}-\beta}{\mbox{erro padrão de }\hat{\beta}}\sim N(0,1)\]

• O intervalo com nível de confiança \((1-\alpha)\times 100\%\) é dado por:

\[\mbox{CI}(\alpha)=\left[\hat{\beta}-z_{\frac{\alpha}{2}}\times\mbox{erro padrão de }\hat{\beta}, \hat{\beta}+z_{\frac{\alpha}{2}}\times\mbox{erro padrão de }\hat{\beta}\right]\]

• Também podemos fazer testes de hipóteses a respeito dos coeficientes, por exemplo, podemos testar se o coeficiente estimado é igual a \(\beta^\star\), na maioria das vezes queremos testar se o coeficiente é igual a zero.

• Com as hipótese que estamos utilizando a estatística de teste para a hipótese nula \(\beta = \beta^\star\) é dada por \(z^\star = \frac{\hat{\beta}-\beta^\star}{\mbox{erro padrão de }\hat{\beta}}\). A distribuição de \(z^\star\) sob a hipótese nula é uma normal padrão.
• Desta forma, podemos calcular o p-valor usando a função de distribuição acumulada da normal padrão.

• Assim como na seção anterior, é comum usar a distribuição \(t\) para fazer testes de hipóteses sobre os coeficientes estimados.
• Com as hipóteses que o termo de erro tem distribuição normal e é homocedastico a estatística \(z^\star\) tem distribuição \(t\) com \(n-2\) graus de liberdade. No caso é \(n-2\) porque estimamos dois parâmetros, \(\alpha\) e \(\beta\).

• Como primeiro exemplo vamos usar o caso das mulheres nos GPs na Índia.

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

fit_women <- lm(water ~ reserved, data= women)

• A função summary() retorna os valores estimados e as estatísticas relevantes.
• A função tidy() do pacote broom retornam essas informações em um tibble.

summary(fit_women)

## 
## Call:
## lm(formula = water ~ reserved, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.991 -14.738  -7.865   2.262 316.009 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   14.738      2.286   6.446 4.22e-10 ***
## reserved       9.252      3.948   2.344   0.0197 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.45 on 320 degrees of freedom
## Multiple R-squared:  0.01688,    Adjusted R-squared:  0.0138 
## F-statistic: 5.493 on 1 and 320 DF,  p-value: 0.0197

library(broom)
tidy(fit_women)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10
## 2 reserved        9.25      3.95      2.34 1.97e- 2

• A função tidy() permite calcular os intervalos de confiança

tidy(fit_women, conf.int=TRUE)

## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10    10.2       19.2
## 2 reserved        9.25      3.95      2.34 1.97e- 2     1.49      17.0

• Intervalo de confiança a 99%

tidy(fit_women, conf.int=TRUE, conf.level=0.99)

## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10    8.81       20.7
## 2 reserved        9.25      3.95      2.34 1.97e- 2   -0.977      19.5

• Repare que esse intervalo de confiança contém o zero.

• Como foi dito, os resultado obtidos com uma variável explicativa valem para modelos com mais variáveis explicativas. Não é objetivo desse curso demonstrar esses resultados em modelos com muitas variáveis (para isso tem a sequência de econometria), mas vamos usar o estudo do salário mínimo como exemplo de regressão múltipla no R.

tidy(fit_minwage, conf.int=TRUE)

## # A tibble: 7 × 7
##   term            estimate std.error statistic p.value conf.low conf.high
##   <chr>              <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 NJ                0.0542    0.0332     1.63  0.103   -0.0111      0.120
## 2 fullPropBefore    0.169     0.0566     2.98  0.00307  0.0574      0.280
## 3 wageBefore        0.0813    0.0389     2.09  0.0374   0.00478     0.158
## 4 chainburgerking  -0.116     0.179     -0.646 0.518   -0.467       0.236
## 5 chainkfc         -0.151     0.183     -0.824 0.411   -0.511       0.209
## 6 chainroys        -0.206     0.187     -1.11  0.270   -0.574       0.161
## 7 chainwendys      -0.220     0.188     -1.17  0.243   -0.591       0.150

• O resultado da regressão é que não podemos rejeitar a hipótese nula que o coeficiente de \(NJ\) é zero, ou seja, o restaurante estar em NJ, onde houve aumento do mínimo, não parece ser relevante para explciar a proporção de empregados em tempo integral.

• Uma das utilidades do modelo de regressão linear é que podemos fazer previsões, de fato, foi assim que apresentamos o modelo na Unidade 4.
• Podemos calcular o erro padrão e construir intervalos de confiança para as previsões do modelo de regressão linear.

• Mais uma vez, por simplicidade, vamos considerar o modelo com uma variável explicativa, \(Y_i = \alpha + \beta X_i + \varepsilon_i\). Estamos interessados em calcular o erro padrão do valor previsto pelo modelo quando \(X\) é igual a um determinado valor \(x\). \[\widehat{Y}=\hat{\alpha}+\hat{\beta}x\]
• Para calcular o variância de \(\widehat{Y}\) devemos ter e mente que \(\hat{\alpha}\) e \(\hat{\beta}\) são possívelmente correlacionados. Quando duas variáveis aleatórias são correlacionadas a variância da soma dessas variáveis não é igual a soma das variãncias.

• Sejam \(X\) e \(Y\) duas variáveis aleatórias. A covariância de \(X\) e \(Y\) é definida por: \[\mbox{Cov}(X,Y)=E[(X-E(X))(Y-E(Y))]=E(XY)-E(X)E(Y)\]
• A correlação, uma versão padronizada da covariância, é dada por: \[\mbox{Cor}(X,Y)=\frac{\mbox{Cov}(X,Y)}{\sqrt{V(X)V(Y)}}\]
• Se duas variáveis aleatórias são independentes uma da outra, então a covariância e a correlação entre elas é zero.

• Para duas variáveis aleatórias, \(X\) e \(Y\), a variância da soma é dada por: \[V(X+Y)=V(X)+V(Y)+2\mbox{Cov}(X,Y)\]
• De forma mais geral: \[V(aX+bY)=a^2V(X)+b^2V(Y)+2ab\mbox{Cov}(X,Y)\]

• Agora podemos calcular a variãncia de \(\widehat{Y}\): \[V(\widehat{Y})=V(\hat{\alpha}+\hat{\beta}x)=V(\hat{\alpha})+V(\hat{\beta})x^2+2x\mbox{Cov}(\hat{\alpha},\hat{\beta})\]

• Erro padrão de \(\widehat{Y}\): \[\mbox{erro padrão de }\widehat{Y}=\sqrt{\widehat{V(\hat{\alpha})}+\widehat{V(\hat{\beta})}x^2+2x\widehat{\mbox{Cov}(\hat{\alpha}, \hat{\beta})}}\]

• Conhecido o erro padrão podemos usar o teorema central do limite para obeter a distribuição do z-score de \(\widehat{Y}\): \[\mbox{z-score de }\widehat{Y}=\frac{\widehat{Y}-(\alpha+\beta x)}{\mbox{erro padrão de }\widehat{Y}} \sim N(0,1)\]

• Para ilustrar esse ponto vamos considerar novamente o exemplo do impacto de ganhar eleições no patrimônimo de um político na Inglaterra.
• O exemplo usou a abrodagem de descontinuidade na regressão, para isso comparou candidatos que ganharam ou perderam por pouco.
• Foram feitas duas regressões para cada partido, uma com os que perderam e outra com os que venceram as eleições.

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

MPs_labour <- filter(MPs, party == "labour")
MPs_tory <- filter(MPs, party == "tory")

labour_fit1 <- lm(ln.net ~ margin, 
                  data = filter(MPs_labour, margin < 0))
labour_fit2 <- lm(ln.net ~ margin, 
                  data = filter(MPs_labour, margin > 0))

tory_fit1 <- lm(ln.net ~ margin, 
                data = filter(MPs_tory, margin < 0))
tory_fit2 <- lm(ln.net ~ margin, 
                data = filter(MPs_tory, margin > 0))

• Como primeiro exercício vamos comparar a riqueza prevista na margem zero em cada regressão para cada partido.

tory_y0 <- augment(tory_fit1, newdata=tibble(margin=0), 
                   interval="confidence", conf.level=0.95)
tory_y1 <- augment(tory_fit2, newdata=tibble(margin=0), 
                   interval="confidence", conf.level=0.95)

tory_y0

## # A tibble: 1 × 4
##   margin .fitted .lower .upper
##    <dbl>   <dbl>  <dbl>  <dbl>
## 1      0    12.5   12.1   13.0

tory_y1

## # A tibble: 1 × 4
##   margin .fitted .lower .upper
##    <dbl>   <dbl>  <dbl>  <dbl>
## 1      0    13.2   12.8   13.6

labour_y0

## # A tibble: 1 × 4
##   margin .fitted .lower .upper
##    <dbl>   <dbl>  <dbl>  <dbl>
## 1      0    12.3   12.0   12.6

labour_y1

## # A tibble: 1 × 4
##   margin .fitted .lower .upper
##    <dbl>   <dbl>  <dbl>  <dbl>
## 1      0    11.9   11.4   12.5

• Repare que no Tory a riqueza prevista quando a margem é zero é maior para os eleitos do que para os não eleitos, o mesmo não acontece com o Labour.
• A questão é saber se essa diferença é significativa. Para os que perderam a riqueza estimada é de 12,5 com intervalo de 95% dado por [12,1, 13]. Para os que ganharam a riqueza estimada é de 13,2 com intervalo de 95% dado por [12,8, 13,6].
• Há uma superposição entre os intervalos.

• Antes de comentar esse resultado, vamos usar um gráfico para ilustrar esse ponto. O primeiro passo é usar a função data_grid() do pacote modelr para criar dois data.frames, um com valores positivos e outro com valores negativos para margem.

y1_range <- data_grid(filter(MPs_tory, margin <=0), margin)
tory_y0 <- augment(tory_fit1, newdata = y1_range, interval="confidence")

y2_range <- data_grid(filter(MPs_tory, margin >= 0), margin)
tory_y1 <- augment(tory_fit2, newdata = y2_range, interval="confidence")

• Com os valores estimados e os intervalos de confiança para margnes positivas e negativas podemos fazer o gráfico. Para destacar a área dos intervalos de confiança usaremos a função geom_ribbon().

ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_ribbon(aes(x=margin, ymin=.lower, ymax=.upper), data = tory_y0, alpha=0.3) +
  geom_line(aes(x=margin, y=.fitted), data = tory_y0) +
  geom_ribbon(aes(x=margin, ymin=.lower, ymax=.upper), data = tory_y1, alpha=0.3) +
  geom_line(aes(x=margin, y=.fitted), data = tory_y1) +
  xlim(-.5, .25) +
  labs(x="margin of victory", y="log net wealth") +
  theme_classic()

• O gráfico ilustra a superposição dos intervalos de confiança.
• Nosso interesse agora é estimar o intervalo de confiança para a diferença entre os valores previstos nos dois modelos para margem igual a zero, pois essa diferença representa o valor estimado para o efeito médio do tratamento em uma abordagem de descontinuidade na regressão.

• Como os modelos foram estimados com dados diferentes, vamos fazer a hipótese que os valores previstos por cada modelo são independentes. Com isso podemos calcular a variância da diferença como a soma das variâncias. \[\mbox{erro padrão de }(\widehat{Y}_1 - \widehat{Y}_0)=\\\sqrt{(\mbox{erro padrão de }\widehat{Y}_1)^2+(\mbox{erro padrão de }\widehat{Y}_0)^2}\]
• O argumento se_fit da função augment() calcula o erro padrão quando definido como TRUE.

tory_y0 <- augment(tory_fit1, newdata=tibble(margin=0), 
                   interval = "confidence", se_fit = TRUE)
tory_y0

## # A tibble: 1 × 5
##   margin .fitted .lower .upper .se.fit
##    <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
## 1      0    12.5   12.1   13.0   0.214

tory_y1 <- augment(tory_fit2, newdata=tibble(margin=0), 
                   interval = "confidence", se_fit = TRUE)
tory_y1

## # A tibble: 1 × 5
##   margin .fitted .lower .upper .se.fit
##    <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
## 1      0    13.2   12.8   13.6   0.192

• Com as informações acima podemo calcular o valor estimado e o erro padrão da diferença.

diff_est <- tory_y1$.fitted - tory_y0$.fitted
diff_est

## [1] 0.6496861

se_diff <- sqrt(tory_y0$.se.fit^2 + tory_y1$.se.fit^2)
se_diff

## [1] 0.2876281

• Intervalo de confiança

CI <- c(diff_est - se_diff * qnorm(0.975), diff_est + se_diff * qnorm(0.975))
CI

## [1] 0.0859455 1.2134268

• Repare que zero não pertence ao intervalo de confiança.

• Teste de hipótese

z_score <- diff_est/se_diff
p.value <- 2*pnorm(abs(z_score), lower.tail = FALSE)
p.value

## [1] 0.02389759

• O p-valor é menor do que 0,05.

• Mesmo com a superposição dos intervalos de confiança a diferença entre as estimativas do dois modelos quando a margem é igual a zero é estatísticamente diferente de zero.
• A análise sugere que ganhar as eleições tem um impacto positivo na riqueza do candidato.