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

• 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 proporçã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 a final de um campeonato). Seja \(p\) a proporção de eleitores que vão votar no candidato de interesse.

• 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 experimentos 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] 1.13755

• 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.0761    1.23      1.31     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] 69 88 89 96 98 86 91 32 53 90 99 56 20 25 95 37 34  5 50 72 40  8 78 59 57
## [26] 80 97  2 17 64 24 22 36 12 82 35 70 28 52  1 30 81 74 19 44 11 45 85 58 49

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.16    1.70    2.86       1
## 2     2 -0.0943  0.00312 0.0974     1
## 3     3  0.503   0.556   0.0524     0
## 4     4 -0.402   3.97    4.38       0
## 5     5  0.330   1.86    1.53       1
## 6     6 -1.12   -0.736   0.383      0

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.16    1.70     2.86        1  1.70   
##  2     2 -0.0943  0.00312  0.0974      1  0.00312
##  3     3  0.503   0.556    0.0524      0  0.503  
##  4     4 -0.402   3.97     4.38        0 -0.402  
##  5     5  0.330   1.86     1.53        1  1.86   
##  6     6 -1.12   -0.736    0.383       0 -1.12   
##  7     7  1.36    1.36     0.00716     0  1.36   
##  8     8  0.0229  1.32     1.30        1  1.32   
##  9     9  0.196  -0.265   -0.460       0  0.196  
## 10    10 -1.31   -0.913    0.399       0 -1.31   
## # ℹ 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.158
## 2     1  1.09

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.158  1.09

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.158    1.09      1.25     0.115

• 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.34756 -0.08601  0.01174  0.01067  0.12137  0.50368

• 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.592257 -0.138294 -0.011056  0.001357  0.147218  0.543378

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

• 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.08 0.211    0.0835

• 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.210   0.201

• 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 processo 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 um candidato 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> 1.0978138, 1.2262773, 1.0781809, 1.4970733, 0.9722847, 1…
## $ se            <dbl> 0.1752122, 0.1962698, 0.2086948, 0.1773209, 0.2037064, 0…
## $ std_error     <dbl> 0.097813810, 0.226277293, 0.078180901, 0.497073283, -0.0…
## $ ci_lower      <dbl> 0.7544042, 0.8415957, 0.6691466, 1.1495307, 0.5730276, 0…
## $ ci_upper      <dbl> 1.441223, 1.610959, 1.487215, 1.844616, 1.371542, 1.6660…
## $ includes_pate <lgl> TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, T…

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

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

pate_sim_coverage(pate_sims_se, level = 0.95)

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

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

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

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

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

• 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.892
## 2   100    0.954
## 3  1000    0.94

• 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 de 3% 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órmula 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.

• Como exemplo vamos analisar os valores previstos e realizados para a proporção de votos em Obama em cada estado nas eleições de 2008.
• Os dados estão em pres08 e polls08, 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

• 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       2
##  3       4
##  4       6
##  5       6
##  6       4
##  7       4
##  8       4
##  9       4
## 10       6

• 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    10 0.02 
## 2       2   110 0.22 
## 3       4   249 0.498
## 4       6   125 0.25 
## 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.02      0.0143  0.00571
## 2       2    0.22      0.229  -0.00857
## 3       4    0.498     0.514  -0.0163 
## 4       6    0.25      0.229   0.0214 
## 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