Lista avaliativa: IC e TH
Probabilidade e Estatística – 2021.2
fnaufel
17/01/2022 (v. 02/02/2022)
- Instruções
- 1 Fumantes
- 2 Graduação em \(4\) anos
- 3 Estimativa do tamanho de uma multidão
- 4 Pulsação média de mulheres
- 5 PlayStations com bootstrapping
- 6 Telefones celulares com bootstrapping
- 7 Sonhando em cores ou em preto e branco
- 8 Professores \(\times\) professoras
- 9 Mulheres \(\times\) maridos
- 10 Temperatura corporal com bootstrapping
- Notas
GABARITO
Instruções
Veja na tabela abaixo o número da questão que foi sorteada para você. Se seu número de matrícula não estiver na tabela, entre em contato comigo pelo Telegram.
matricula
questao
119060007
1
119060025
3
119060038
8
119060043
9
120060010
4
120060014
6
120060029
7
120060034
10
216060055
5
219060053
2
Clique o botão
Code, no início desta página, para baixar o arquivo Rmd deste documento.Certifique-se de que o pacote
fnaufelRmdestá instalado. Use os comandosinstall.packages("devtools") devtools::install_github("fnaufel/fnaufelRmd")Edite o arquivo Rmd para resolver a sua questão.
Cada item que pede para você construir um intervalo de confiança ou para você fazer um teste de hipóteses deve ser resolvido de duas maneiras (exceto quando houver instrução dizendo o contrário):
Fazendo os cálculos passo a passo, como nos vídeos do curso, e
Usando as funções prontas do R para IC e testes de hipótese que foram usadas nos vídeos do curso.
Escreva o máximo possível sobre o seu raciocínio. Justifique suas respostas.
Teste suas respostas. Tudo deve estar executando sem erros.
Envie no Moodle:
O arquivo Rmd com as suas resoluções e
Um vídeo de até \(5\) minutos explicando as suas resoluções.
Bom trabalho.
1 Fumantes
Uma amostra aleatória de \(86\) pessoas inclui \(41\) fumantes.
Construa um intervalo de confiança de \(95\%\) para o percentual de fumantes na população de onde foi coletada a amostra.
Condições
A distribuição da população é normal e \(n \geq 10\)?
Não sabemos se a distribuição da população é normal.
Caso a distribuição da população não seja conhecida, \(n \geq 30\)?
Sim.
A amostra é aleatória, representativa da população?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
O tamanho da amostra é menos do que \(10\%\) da população?
Sim.
Há pelo menos 10 fracassos e pelo menos 10 sucessos na amostra?
Sim.
Passo a passo
n <- 86 sucessos <- 41 nivel_confianca <- .95 pchapeu <- sucessos / n cat('pchapeu =', pchapeu) ## pchapeu = 0,4767442 erro_padrao <- sqrt(pchapeu * (1 - pchapeu) / n) cat('\nEP =', erro_padrao) ## ## EP = 0,05385804 valor_critico <- -qnorm((1 - nivel_confianca)/2) cat('\nValor crítico =', valor_critico) ## ## Valor crítico = 1,959964 margem_erro <- erro_padrao * valor_critico cat('\nME =', margem_erro) ## ## ME = 0,1055598 intervalo <- pchapeu + c(-1, 1) * margem_erro cat('\nIC = [', intervalo[1], ';', intervalo[2], ']') ## ## IC = [ 0,3711844 ; 0,582304 ]Graficamente, segundo esta amostra, a distribuição amostral de \(\hat p\) é a normal abaixo, com média \(0{,}477\) e desvio padrão \(0{,}054\). A área preenchida corresponde a \(95\%\) de probabilidade:
Faça um teste de hipóteses para avaliar a afirmação de que o percentual de fumantes na população é diferente de \(50\%\). Use \(\alpha = 0{,}05\).
Passo a passo
n <- 86 sucessos <- 41 # Valor de p na hipótese nula pzero <- .5 # Supondo a hipótese nula erro_padrao <- sqrt(pzero * (1 - pzero) / n) cat('\nEP =', erro_padrao) ## ## EP = 0,05391639 pchapeu <- sucessos / n cat('\npchapeu =', pchapeu) ## ## pchapeu = 0,4767442 # Valor p valor_p <- 2 * pnorm( pchapeu, mean = pzero, sd = erro_padrao ) cat('\nValor p =', valor_p) ## ## Valor p = 0,6662276Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição normal, com média \(0{,}500\) e desvio padrão \(\approx 0{,}054\):
Note que esta é a distribuição amostral de \(\hat p\) supondo que \(H_0\) é verdadeira; sua média é \(0{,}5\), diferente da distribuição amostral de \(\hat p\) baseada na amostra que temos, desenhada acima, cuja média era \(0{,}477\).
Como o valor \(p\) \((\approx 0{,}666)\) é muito maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Usando a função
prop.test()n <- 86 sucessos <- 41 prop.test(sucessos, n, p = .5, correct = FALSE) ## ## 1-sample proportions test without continuity correction ## ## data: sucessos out of n, null probability 0.5 ## X-squared = 0,18605, df = 1, p-value = 0,6662 ## alternative hypothesis: true p is not equal to 0,5 ## 95 percent confidence interval: ## 0,3744554 0,5810217 ## sample estimates: ## p ## 0,4767442
2 Graduação em \(4\) anos
Qual o tamanho mínimo da amostra necessária para estimar, com um nível de confiança de \(95\%\) e uma margem de erro de \(0{,}05\) (\(5\) pontos percentuais), a proporção de alunos universitários que completam a graduação em \(4\) anos ou menos? Nada mais é sabido sobre a população.
Explique por que você escolheu o valor de \(\hat p\) que você escolheu.
Quando construímos um intervalo de confiança com nível de confiança de \(95\%\), a margem de erro é calculada como
\[ \text{ME} = z^\star \cdot \sqrt{\frac{\hat p (1 - \hat p)}{n}} \]
onde \(z^\star = {}\)
-qnorm((1 - .95)/2)\({} \approx 1{,}96\).Se isolarmos \(n\), chegamos a
\[ n = \hat p (1-\hat p) \left(\frac{z^\star}{\text{ME}}\right)^2 \]
Não sabemos o valor de \(\hat p\). Vamos pensar no maior \(n\) possível, que corresponde a \(\hat p = 0{,}5\). Então,
\[ n = \hat p (1-\hat p) \left(\frac{z^\star}{\text{ME}}\right)^2 = 0{,}5 \cdot 0{,}5 \cdot \left( \frac{1{,}96}{0{,}05} \right)^2 \approx 385 \]
Ou seja, na pior das hipóteses, vamos precisar de uma amostra de \(385\) ou mais pessoas.
Em uma amostra de \(400\) alunos universitários, \(150\) completaram a graduação em \(4\) anos ou menos. Teste a afirmação de que, na população, menos de \(40\%\) dos alunos universitários completam a graduação neste tempo.
Condições
A distribuição da população é normal e \(n \geq 10\)?
Não sabemos se a distribuição da população é normal.
Caso a distribuição da população não seja conhecida, \(n \geq 30\)?
Sim.
A amostra é aleatória, representativa da população?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
O tamanho da amostra é menos do que \(10\%\) da população?
Sim.
Há pelo menos 10 fracassos e pelo menos 10 sucessos na amostra?
Sim.
Passo a passo
n <- 400 sucessos <- 150 # Valor de p na hipótese nula pzero <- .4 # Supondo a hipótese nula erro_padrao <- sqrt(pzero * (1 - pzero) / n) cat('\nEP =', erro_padrao) ## ## EP = 0,0244949 pchapeu <- sucessos / n cat('\npchapeu =', pchapeu) ## ## pchapeu = 0,375 # Valor p valor_p <- pnorm( pchapeu, mean = pzero, sd = erro_padrao ) cat('\nValor p =', valor_p) ## ## Valor p = 0,1537171Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição normal, com média \(0{,}400\) e desvio padrão \(\approx 0{,}024\):
Como o valor \(p\) \((\approx 0{,}154)\) é maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Como a hipótese alternativa é \(H_A : p < 0{,}4\), o teste é unilateral.
Usando a função
prop.test()n <- 400 sucessos <- 150 prop.test(sucessos, n, p = .4, alternative = 'less', correct = FALSE) ## ## 1-sample proportions test without continuity correction ## ## data: sucessos out of n, null probability 0.4 ## X-squared = 1,0417, df = 1, p-value = 0,1537 ## alternative hypothesis: true p is less than 0,4 ## 95 percent confidence interval: ## 0,0000000 0,4155303 ## sample estimates: ## p ## 0,375
3 Estimativa do tamanho de uma multidão1
Observe a foto abaixo, que mostra parte de uma multidão. A quantidade de pessoas em cada setor quadrado foi estimada por um algoritmo.
Os valores obtidos para os setores da foto foram
setores <- c(
200, 400, 800, 1000, 400, 800, 700, 1000, 400, 200, 700, 900, 400, 400, 700, 900,
600, 700, 900, 650, 672, 700, 900, 450, 566, 700, 900
)Construa um intervalo de confiança de \(95\%\) para a quantidade média de pessoas por setor na multidão inteira (inclusive as partes que não aparecem na foto).
A população consiste de todas as quantidades de pessoas por setor na multidão inteira.
O vetor
setoresé uma amostra.Como não sabemos o desvio padrão da população, vamos usar o teste \(t\) para estimar a média de pessoas por setor.
Condições
A amostra é aleatória, representativa da população?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
A amostra vem de uma população aproximadamente normal?
Podemos fazer um teste de normalidade:
setores %>% shapiro.test() ## ## Shapiro-Wilk normality test ## ## data: . ## W = 0,93376, p-value = 0,08547O valor \(p\) acima de \(0{,}05\) é evidência de que os dados seguem uma distribuição normal.
Podemos fazer um histograma:
tibble(n = setores) %>% ggplot(aes(n)) + geom_histogram(bins = 5) + labs(y = NULL)
A forma do histograma é aproximadamente normal.
Passo a passo
n = length(setores) cat('n =', n) ## n = 27 xbarra <- mean(setores) cat('\nMédia amostral =', xbarra) ## ## Média amostral = 653,2593 s <- sd(setores) cat('\nDesvio padrão amostral =', s) ## ## Desvio padrão amostral = 230,7807 erro_padrao <- s / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 44,41377 gl <- n - 1 cat('\nGraus de liberdade = ', gl) ## ## Graus de liberdade = 26 nivel_confianca <- .95 valor_critico <- -qt((1 - nivel_confianca)/2, df = gl) cat('\nValor crítico =', valor_critico) ## ## Valor crítico = 2,055529 margem_erro <- erro_padrao * valor_critico cat('\nME =', margem_erro) ## ## ME = 91,29381 intervalo <- xbarra + c(-1, 1) * margem_erro cat('\nIC = [', intervalo[1], ';', intervalo[2], ']') ## ## IC = [ 561,9654 ; 744,5531 ]Graficamente, o intervalo de confiança corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(26\) graus de liberdade:
Observe que esta distribuição \(t\) corresponde à variável padronizada
\[ \frac{x - \bar x}{s/\sqrt{n}} \]
Trabalhamos com esta variável padronizada porque a distribuição \(t\) sempre tem média zero.
Aqui, cada valor no eixo \(x\) corresponde ao número de desvios padrão acima ou abaixo da média.
De fato, o limite inferior da área preenchida corresponde a
(intervalo[1] - xbarra) / erro_padrao ## [1] -2,055529e o limite superior da área preenchida corresponde a
(intervalo[2] - xbarra) / erro_padrao ## [1] 2,055529Usando o resultado do item anterior e considerando que a multidão inteira cobre \(50\) setores, construa um intervalo de confiança para o total de pessoas na multidão.
Basta multiplicar os valores do intervalo de confiança por \(50\):
intervalo * 50 ## [1] 28098,27 37227,65Faça um teste de hipóteses para avaliar a afirmação de que a multidão inteira consiste de mais de \(35\) mil pessoas. Use \(\alpha = 0{,}05\).
A multidão inteira ter \(35\) mil pessoas equivale a
\[ \mu = \frac{35.000}{50} = 700 \]
onde \(\mu\) é a média de pessoas por setor.
Antes de fazer qualquer cálculo, percebemos que o valor \(700\) está dentro do intervalo de confiança de \(95\%\) para a média de pessoas por setor.
Isto significa que a hipótese de nulidade não vai ser rejeitada.
Passo a passo
n <- length(setores) cat('n =', n) ## n = 27 s <- sd(setores) cat('\nDesvio padrão amostral =', s) ## ## Desvio padrão amostral = 230,7807 xbarra <- mean(setores) cat('\nMédia amostral =', xbarra) ## ## Média amostral = 653,2593 # Valor de μ na hipótese nula mizero <- 700 # Supondo a hipótese nula erro_padrao <- s / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 44,41377 # Média amostral padronizada (supondo a hipótese nula) tzero <- (xbarra - mizero) / erro_padrao cat('\ntzero =', tzero) ## ## tzero = -1,052393 # Graus de liberdade gl <- n - 1 # Valor p valor_p <- pt( tzero, df = gl, lower.tail = FALSE ) cat('\nValor p =', valor_p) ## ## Valor p = 0,8488489Como o valor \(p\) \((\approx 0{,}849)\) é muito maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(26\) graus de liberdade, que corresponde à variável padronizada
\[ \frac{x - \bar x}{s/\sqrt{n}} \]
Como a hipótese alternativa é \(H_A : \mu > 700\), o teste é unilateral.
Usando a função
t.test()t.test(x = setores, alternative = 'greater', mu = mizero) ## ## One Sample t-test ## ## data: setores ## t = -1,0524, df = 26, p-value = 0,8488 ## alternative hypothesis: true mean is greater than 700 ## 95 percent confidence interval: ## 577,5063 Inf ## sample estimates: ## mean of x ## 653,2593
4 Pulsação média de mulheres
Qual o tamanho mínimo da amostra necessária para estimar, com um nível de confiança de \(95\%\) e uma margem de erro de \(2\) batimentos por minuto (bpm), a pulsação média de mulheres adultas? De estudos anteriores, podemos usar a informação de que o desvio padrão da população é \(\sigma = 12{,}54\) bpm.
Quando construímos um intervalo de confiança com nível de confiança de \(95\%\), a margem de erro é calculada como
\[ \text{ME} = z^\star \cdot \frac{\sigma}{\sqrt{n}} \]
onde \(z^\star = {}\)
-qnorm((1 - .95)/2)\({} \approx 1{,}96\).Se isolarmos \(n\), chegamos a
\[ n = \left(\frac{z^\star \cdot \sigma}{\text{ME}}\right)^2 \]
Substituindo os valores que temos:
\[ n = \left(\frac{1{,}96 \cdot 12{,}54}{2}\right)^2 \approx 152 \]
Vamos precisar de uma amostra de \(152\) ou mais pessoas.
O vetor
pulsacaocontém os valores das pulsações (em bpm) de diversas mulheres.pulsacao <- c( 80, 94, 58, 66, 56, 82, 78, 86, 88, 56, 36, 66, 84, 76, 78, 64, 66, 78, 60, 64, 84, 82, 70, 74, 86, 90, 88, 90, 90, 94, 68, 90, 82, 80, 74, 56, 100, 74, 76, 76, 86, 74, 66, 62, 78, 68, 62, 62, 74, 104, 54, 74, 74, 84, 60, 52, 84, 66, 56, 66, 84, 64, 64, 78, 104, 84, 84, 62, 82, 64, 78, 76, 86, 72, 82, 82, 66, 58, 62, 80, 58, 82, 76, 70, 88, 44, 58, 64, 68, 70, 58, 72, 80, 70, 74, 72, 66, 82, 98, 74, 80, 82, 56, 78, 100, 56, 68, 72, 64, 94, 72, 56, 80, 76, 56, 56, 54, 72, 70, 80, 72, 62, 96, 72, 60, 98, 72, 80, 84, 68, 86, 72, 80, 64, 82, 80, 62, 74, 96, 72, 96, 64, 62, 82, 76, 94, 74 )Teste a afirmação de que a pulsação média de mulheres adultas é menor do que \(70\) bpm. Use \(\alpha = 0{,}05\). Use o fato de que \(\sigma = 12{,}54\) bpm.
Condições
Como temos o valor do desvio padrão da população, podemos usar o teste \(z\), cujas condições são:
A amostra é aleatória, representativa da população?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
O tamanho da amostra é menor do que \(10\%\) da população?
A população é o conjunto de todas as mulheres adultas (da cidade? do país?)
Nossa amostra tem \(147\), que é menos do que \(10\%\) da população.
Passo a passo
n <- length(pulsacao) cat('n =', n) ## n = 147 xbarra <- mean(pulsacao) cat('\nMédia amostral =', xbarra) ## ## Média amostral = 74,04082 # Desvio padrão da população (dado) sigma <- 12.54 # Valor de μ na hipótese nula mizero <- 70 erro_padrao <- sigma / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 1,034282 # Valor p valor_p <- pnorm( xbarra, mean = mizero, sd = erro_padrao ) cat('\nValor p =', valor_p) ## ## Valor p = 0,9999533Como o valor \(p\) \((\approx 1{,}000)\) é muito maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Na verdade, como a média da amostra é maior do que \(70\)bpm, esta amostra não ajuda a refutar a hipótese de nulidade, que diz que \(\mu < 70\).
Além disso, o valor baixo do erro padrão diz que as amostras não variam muito em suas médias; ou seja, esta é uma amostra típica.
Por isso, é praticamente impossível refutar a hipótese de nulidade com esta amostra.
Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da normal na seguinte figura. Perceba como a área à esquerda do valor de \(\bar x = 74{,}041\) é praticamente \(100\%\) da área inteira abaixo da curva:
Como a hipótese alternativa é \(H_A : \mu < 70\), o teste é unilateral.
Usando a função
z.test()library(BSDA) ## Loading required package: lattice z.test( pulsacao, alternative = 'less', mu = 70, sigma.x = sigma ) ## ## One-sample z-Test ## ## data: pulsacao ## z = 3,9069, p-value = 1 ## alternative hypothesis: true mean is less than 70 ## 95 percent confidence interval: ## NA 75,74206 ## sample estimates: ## mean of x ## 74,04082Refaça o item anterior sem usar o fato de que \(\sigma = 12{,}54\) bpm.
Condições
Como agora não temos o desvio padrão da população, precisamos fazer um teste \(t\), cujas condições são:
A amostra é aleatória, representativa da população?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
A amostra vem de uma população aproximadamente normal?
Podemos fazer um teste de normalidade:
pulsacao %>% shapiro.test() ## ## Shapiro-Wilk normality test ## ## data: . ## W = 0,98912, p-value = 0,3102O valor \(p\) acima de \(0{,}05\) é evidência de que os dados seguem uma distribuição normal.
Podemos fazer um histograma:
tibble(n = pulsacao) %>% ggplot(aes(n)) + geom_histogram(bins = 9) + labs(y = NULL)
A forma do histograma é aproximadamente normal.
Passo a passo
n <- length(pulsacao) cat('n =', n) ## n = 147 s <- sd(pulsacao) cat('\nDesvio padrão amostral =', s) ## ## Desvio padrão amostral = 12,54356 xbarra <- mean(pulsacao) cat('\nMédia amostral =', xbarra) ## ## Média amostral = 74,04082 # Valor de μ na hipótese nula mizero <- 70 # Supondo a hipótese nula erro_padrao <- s / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 1,034575 # Média amostral padronizada (supondo a hipótese nula) tzero <- (xbarra - mizero) / erro_padrao cat('\ntzero =', tzero) ## ## tzero = 3,905774 # Graus de liberdade gl <- n - 1 cat('\ngl =', gl) ## ## gl = 146 # Valor p valor_p <- pt( tzero, df = gl ) cat('\nValor p =', valor_p) ## ## Valor p = 0,9999284Como o valor \(p\) \((\approx 1{,}000)\) é muito maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(146\) graus de liberdade, que corresponde à variável padronizada
\[ \frac{x - \bar x}{s/\sqrt{n}} \]
Trabalhamos com esta variável padronizada porque a distribuição \(t\) sempre tem média zero.
Aqui, cada valor no eixo \(x\) corresponde ao número de desvios padrão acima ou abaixo da média.
Como a hipótese alternativa é \(H_A : \mu < 70\), o teste é unilateral.
Usando a função
t.test()t.test(x = pulsacao, alternative = 'less', mu = mizero) ## ## One Sample t-test ## ## data: pulsacao ## t = 3,9058, df = 146, p-value = 0,9999 ## alternative hypothesis: true mean is less than 70 ## 95 percent confidence interval: ## -Inf 75,75341 ## sample estimates: ## mean of x ## 74,04082Compare e comente os resultados.
Observe que
O desvio padrão populacional \(\sigma = 12{,}540\) é bem próximo do desvio padrão amostral \(s = 12{,}544\).
Com \(146\) graus de liberdade, a distribuição \(t\) é muito próxima da normal padrão:
Estas semelhanças tornam os cálculos dos valores \(p\) bem parecidos, tanto no teste \(z\) quanto no teste \(t\).
5 PlayStations com bootstrapping
Em uma amostra \(A\) de \(200\) pessoas, \(124\) possuem um PlayStation.
O método de bootstrapping consiste em simular muitas amostras — com reposição — a partir dos elementos desta amostra original \(A\). A partir destas muitas amostras simuladas, você pode calcular intervalos de confiança como explicado abaixo:
Armazene, em uma lista, \(10\) mil amostras, cada uma com \(200\) pessoas sorteadas com reposição a partir da amostra original \(A\) acima (com todas as \(200\) pessoas tendo a mesma probabilidade de sorteio).
Poderíamos criar um vetor com os \(200\) elementos da amostra original para passar para a função
sample:amostra_original <- rep(c(TRUE, FALSE), c(124, 76))Mas podemos pensar de outro modo: ao escolhermos um elemento da amostra original, a probabilidade de obtermos
TRUEé igual a \(124/200\).Ou seja, sortear \(200\) pessoas desta amostra com reposição equivale a sortear \(200\) valores booleanos, com
TRUEtendo probabilidade \(124/200\) de ser sorteado a cada vez:n_amostras <- 1e4 n <- 200 sucessos <- 124 p <- sucessos / n sortear_amostra_bool <- function(n, p) { sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(p, 1 - p)) }Para construir a lista de \(10\) mil amostras:
amostras <- replicate( n_amostras, sortear_amostra_bool(n, p), simplify = FALSE )Vamos examinar a lista gerada:
str(amostras) ## List of 10000 ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] FALSE FALSE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE FALSE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] FALSE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE FALSE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE FALSE FALSE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE FALSE ... ## $ : logi [1:200] FALSE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE FALSE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] FALSE TRUE FALSE FALSE ... ## $ : logi [1:200] FALSE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] FALSE TRUE FALSE FALSE ... ## $ : logi [1:200] FALSE FALSE FALSE FALSE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] FALSE FALSE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE FALSE ... ## $ : logi [1:200] FALSE TRUE FALSE FALSE ... ## $ : logi [1:200] FALSE TRUE TRUE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] FALSE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE TRUE FALSE ... ## $ : logi [1:200] TRUE FALSE TRUE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE TRUE FALSE ... ## $ : logi [1:200] TRUE TRUE FALSE TRUE ... ## $ : logi [1:200] TRUE FALSE FALSE TRUE ... ## [list output truncated]Armazene, em um vetor de \(10\) mil elementos, as proporções de possuidores de PlayStations nas amostras obtidas no item anterior.
Lembre-se de que calcular a proporção de valores
TRUEem um vetor booleano equivale a calcular a média (já queTRUEé visto como \(1\), eFALSEé visto como \(0\)).As funções da família
mapaplicam uma operação (aqui, a média) a cada elemento de um vetor, retornando uma lista. Veja?purrr::mappara mais detalhes sobre programação funcional em R.proporcoes <- amostras %>% map_dbl(~mean(.))str(proporcoes) ## num [1:10000] 0,615 0,68 0,59 0,615 0,635 0,61 0,575 0,63 ...Imprima os valores deste vetor que estão no quantil \(0{,}05\) e no quantil \(0{,}95\). Estes são os extremos do intervalo de confiança para a proporção de possuidores de PlayStations na população, com \(95\%\) de confiança, calculados via bootstrapping.
Uma das vantagens do bootstrapping é não precisar de qualquer conhecimento sobre a população — sua distribuição, sua variância etc.
ic_bs <- proporcoes %>% quantile(c(.05, .95)) ic_bs ## 5% 95% ## 0,565 0,675Use a função
prop.test()do R para construir um intervalo de confiança para a proporção de possuidores de PlayStations na população, com \(95\%\) de confiança — não é necessário fazer os cálculos passo a passo.ic_teste <- prop.test(sucessos, n, correct = FALSE)$conf.int ic_teste ## [1] 0,5510672 0,6844099 ## attr(,"conf.level") ## [1] 0,95Compare com o resultado do bootstrapping.
Os centros dos intervalos são
Bootstrap: \(0{,}620\)
Teste: \(0{,}618\)
Os extremos de um intervalo são próximos dos extremos do outro intervalo.
Como o bootstrapping consiste em sortear elementos da amostra original, os resultados podem ser diferentes a cada execução, mas os resultados costumam ser próximos aos dos testes paramétricos — como
prop.test.Se tomarmos uma quantidade maior de amostras — aqui, foram \(10\) mil — a precisão melhorará mais ainda. Experimente.
6 Telefones celulares com bootstrapping
Uma amostra das quantidades de radiação emitidas por telefones celulares consiste dos seguintes valores:
radiacao <- c(38, 55, 86, 91, 99, 103, 145)O método de bootstrapping consiste em simular muitas amostras — com reposição — a partir dos elementos desta amostra original. A partir destas muitas amostras simuladas, você pode calcular intervalos de confiança como explicado abaixo:
Armazene, em uma lista, \(10\) mil amostras, cada uma com \(7\) elementos sorteados com reposição a partir do vetor
radiacao(com todos os elementos tendo a mesma probabilidade de sorteio).n_amostras <- 1e4 n <- length(radiacao) amostras <- replicate( n_amostras, sample(radiacao, n, replace = TRUE), simplify = FALSE )Vamos examinar a lista gerada:
str(amostras) ## List of 10000 ## $ : num [1:7] 91 103 86 99 86 99 55 ## $ : num [1:7] 145 55 99 38 86 38 86 ## $ : num [1:7] 86 86 99 55 99 145 38 ## $ : num [1:7] 103 55 86 38 86 145 38 ## $ : num [1:7] 55 145 86 86 99 91 145 ## $ : num [1:7] 145 55 103 145 55 145 145 ## $ : num [1:7] 86 99 55 99 99 86 99 ## $ : num [1:7] 99 55 99 91 86 91 103 ## $ : num [1:7] 91 38 99 86 103 99 145 ## $ : num [1:7] 86 55 91 38 91 91 86 ## $ : num [1:7] 38 145 99 55 103 99 55 ## $ : num [1:7] 99 86 145 55 103 99 91 ## $ : num [1:7] 38 103 99 86 99 38 145 ## $ : num [1:7] 38 55 103 86 86 86 99 ## $ : num [1:7] 38 86 91 103 38 86 103 ## $ : num [1:7] 145 38 103 38 145 103 38 ## $ : num [1:7] 86 91 103 55 103 86 55 ## $ : num [1:7] 55 55 55 91 103 86 38 ## $ : num [1:7] 86 86 38 145 86 38 86 ## $ : num [1:7] 99 99 145 91 86 55 103 ## $ : num [1:7] 86 38 91 91 38 145 86 ## $ : num [1:7] 99 38 55 38 91 103 103 ## $ : num [1:7] 38 145 55 91 103 38 99 ## $ : num [1:7] 91 38 86 91 99 55 38 ## $ : num [1:7] 38 38 145 55 86 145 91 ## $ : num [1:7] 145 91 145 103 103 91 55 ## $ : num [1:7] 86 103 86 99 86 99 55 ## $ : num [1:7] 145 38 103 55 86 38 99 ## $ : num [1:7] 103 103 86 86 103 91 91 ## $ : num [1:7] 103 103 38 86 91 103 38 ## $ : num [1:7] 91 38 38 38 145 145 99 ## $ : num [1:7] 86 55 55 38 99 145 91 ## $ : num [1:7] 86 55 91 91 91 86 38 ## $ : num [1:7] 103 86 86 91 145 103 86 ## $ : num [1:7] 145 55 145 103 55 86 99 ## $ : num [1:7] 86 91 38 145 55 145 38 ## $ : num [1:7] 99 55 99 99 86 103 99 ## $ : num [1:7] 55 145 103 103 86 99 86 ## $ : num [1:7] 145 145 145 99 91 99 99 ## $ : num [1:7] 91 145 55 145 145 91 91 ## $ : num [1:7] 103 99 38 91 38 86 91 ## $ : num [1:7] 145 91 103 86 38 55 145 ## $ : num [1:7] 91 38 55 103 99 38 91 ## $ : num [1:7] 55 38 103 38 145 55 145 ## $ : num [1:7] 91 38 103 103 86 55 55 ## $ : num [1:7] 55 99 145 99 38 103 91 ## $ : num [1:7] 38 103 86 103 38 99 145 ## $ : num [1:7] 103 145 38 91 86 38 103 ## $ : num [1:7] 38 86 99 86 38 103 86 ## $ : num [1:7] 145 103 91 91 99 99 38 ## $ : num [1:7] 86 103 145 86 86 86 55 ## $ : num [1:7] 91 103 55 99 99 103 91 ## $ : num [1:7] 86 145 91 55 145 38 38 ## $ : num [1:7] 99 91 86 91 145 99 145 ## $ : num [1:7] 55 103 99 99 145 86 145 ## $ : num [1:7] 99 103 91 55 55 91 103 ## $ : num [1:7] 145 91 55 38 91 55 55 ## $ : num [1:7] 38 99 99 86 38 91 99 ## $ : num [1:7] 91 38 38 91 55 55 86 ## $ : num [1:7] 55 145 86 99 99 86 91 ## $ : num [1:7] 86 38 86 86 55 99 103 ## $ : num [1:7] 86 145 91 91 91 99 55 ## $ : num [1:7] 103 38 86 145 38 91 91 ## $ : num [1:7] 86 145 55 86 103 55 99 ## $ : num [1:7] 99 86 55 38 38 38 99 ## $ : num [1:7] 145 91 91 91 99 86 103 ## $ : num [1:7] 55 145 91 38 145 145 99 ## $ : num [1:7] 55 91 38 38 38 55 103 ## $ : num [1:7] 38 145 99 145 145 99 55 ## $ : num [1:7] 55 55 55 86 103 55 86 ## $ : num [1:7] 103 145 145 86 91 145 91 ## $ : num [1:7] 99 86 55 103 38 86 55 ## $ : num [1:7] 145 99 38 145 55 38 86 ## $ : num [1:7] 55 91 91 99 99 91 86 ## $ : num [1:7] 55 145 86 99 145 103 86 ## $ : num [1:7] 86 99 91 86 145 103 55 ## $ : num [1:7] 86 145 103 145 103 38 86 ## $ : num [1:7] 145 99 55 145 103 103 86 ## $ : num [1:7] 103 55 99 145 103 103 86 ## $ : num [1:7] 55 103 91 55 99 38 145 ## $ : num [1:7] 99 145 38 55 99 145 38 ## $ : num [1:7] 91 86 91 55 86 99 145 ## $ : num [1:7] 38 55 91 55 145 145 55 ## $ : num [1:7] 55 55 91 91 99 55 99 ## $ : num [1:7] 55 145 99 91 145 103 103 ## $ : num [1:7] 55 99 145 91 38 55 91 ## $ : num [1:7] 99 145 99 38 145 91 99 ## $ : num [1:7] 103 55 145 55 86 103 103 ## $ : num [1:7] 103 99 91 86 99 99 86 ## $ : num [1:7] 55 38 55 38 103 86 86 ## $ : num [1:7] 99 91 38 55 55 91 86 ## $ : num [1:7] 103 145 145 38 91 99 91 ## $ : num [1:7] 99 38 103 91 145 38 103 ## $ : num [1:7] 91 145 55 91 55 91 55 ## $ : num [1:7] 103 145 55 38 55 91 103 ## $ : num [1:7] 91 86 145 55 145 91 99 ## $ : num [1:7] 103 91 91 86 145 86 38 ## $ : num [1:7] 86 38 86 55 103 91 38 ## $ : num [1:7] 38 91 145 91 38 86 86 ## [list output truncated]Armazene, em um vetor de \(10\) mil elementos, as médias das amostras obtidas no item anterior.
As funções da família
mapaplicam uma operação (aqui, a média) a cada elemento de um vetor, retornando uma lista. Veja?purrr::mappara mais detalhes sobre programação funcional em R.medias <- amostras %>% map_dbl(~mean(.))str(medias) ## num [1:10000] 88,4 78,1 86,9 78,7 ...Imprima os valores deste vetor que estão no quantil \(0{,}05\) e no quantil \(0{,}95\). Estes são os extremos do intervalo de confiança para a radiação média na população, com \(95\%\) de confiança, calculados via bootstrapping.
Uma das vantagens do bootstrapping é não precisar de qualquer conhecimento sobre a população — sua distribuição, sua variância etc.
ic_bs <- medias %>% quantile(c(.05, .95)) ic_bs ## 5% 95% ## 68,57143 108,57143Use a função
t.test()do R para construir um intervalo de confiança para a radiação média, com \(95\%\) de confiança — não é necessário fazer os cálculos passo a passo.ic_teste <- t.test(radiacao)$conf.int ic_teste ## [1] 56,11225 120,17346 ## attr(,"conf.level") ## [1] 0,95Compare com o resultado do bootstrapping.
Os centros dos intervalos são
Bootstrap: \(88{,}571\)
Teste: \(88{,}143\)
Os extremos de um intervalo são mais ou menos próximos dos extremos do outro intervalo. Neste exemplo, a amostra original era pequena; se fosse maior, os intervalos seriam mais parecidos.
Como o bootstrapping consiste em sortear elementos da amostra original, os resultados podem ser diferentes a cada execução, mas os resultados costumam ser próximos aos dos testes paramétricos — como
prop.test.
7 Sonhando em cores ou em preto e branco2
Dentre \(306\) pessoas com mais de \(55\) anos, \(68\) disseram que sonham em preto e branco.
Dentre \(298\) pessoas com menos de \(25\) anos, \(13\) disseram que sonham em preto e branco.
Construa um intervalo de confiança de \(99\%\) para a diferença entre as duas proporções.
Condições
A distribuição da população é normal e \(n \geq 10\)?
Não sabemos se a distribuição da população é normal.
Caso a distribuição da população não seja conhecida, \(n \geq 30\)?
Sim.
As amostras são aleatórias, representativas das populações?
Supomos que sim.
Os valores da amostra são independentes entre si?
Supomos que sim.
O tamanho da amostra é menos do que \(10\%\) da população?
Sim.
Há pelo menos 10 fracassos e pelo menos 10 sucessos em cada amostra?
Sim.
Os grupos são independentes?
Supomos que sim (não são parentes nem de outra forma relacionados).
Passo a passo
n1 <- 306 sucessos1 <- 68 pchapeu1 <- sucessos1 / n1 cat('\npchapeu1 =', pchapeu1) ## ## pchapeu1 = 0,2222222 n2 <- 298 sucessos2 <- 13 pchapeu2 <- sucessos2 / n2 cat('\npchapeu2 =', pchapeu2) ## ## pchapeu2 = 0,04362416 dchapeu <- pchapeu1 - pchapeu2 cat('\ndchapeu =', dchapeu) ## ## dchapeu = 0,1785981 nivel_confianca <- .99 erro_padrao <- sqrt( (pchapeu1 * (1 - pchapeu1) / n1) + (pchapeu2 * (1 - pchapeu2) / n2) ) cat('\nEP =', erro_padrao) ## ## EP = 0,0265488 valor_critico <- -qnorm((1 - nivel_confianca)/2) cat('\nValor crítico =', valor_critico) ## ## Valor crítico = 2,575829 margem_erro <- erro_padrao * valor_critico cat('\nME =', margem_erro) ## ## ME = 0,06838517 intervalo <- dchapeu + c(-1, 1) * margem_erro cat('\nIC = [', intervalo[1], ';', intervalo[2], ']') ## ## IC = [ 0,1102129 ; 0,2469832 ]Graficamente, segundo esta amostra, a distribuição amostral de \(\hat d\) (a diferença estimada entre as proporções) é a normal abaixo, com média \(0{,}179\) e desvio padrão \(0{,}027\). A área preenchida corresponde a \(99\%\) de probabilidade:
Usando a função
prop.testprop.test( matrix( c( sucessos1, sucessos2, n1 - sucessos1, n2 - sucessos2 ), ncol = 2 ), conf.level = .99, correct = FALSE ) ## ## 2-sample test for equality of proportions without continuity correction ## ## data: matrix(c(sucessos1, sucessos2, n1 - sucessos1, n2 - sucessos2), ncol = 2) ## X-squared = 41,471, df = 1, p-value = 0,0000000001197 ## alternative hypothesis: two.sided ## 99 percent confidence interval: ## 0,1102129 0,2469832 ## sample estimates: ## prop 1 prop 2 ## 0,22222222 0,04362416Teste a afirmação de que a proporção de pessoas com mais de \(55\) anos que sonham em preto e branco é maior do que a proporção de pessoas com menos de \(25\) anos que sonham em preto e branco. Use \(\alpha = 0{,}01\).
Como o intervalo de confiança para a diferença entre as proporções, calculado no item anterior, não inclui o zero, e como \(\hat p_1 > \hat p _2\), já sabemos, antes de qualquer cálculo, que a hipótese de nulidade vai ser rejeitada.
Neste nível de confiança, a amostra é evidência de que a proporção de pessoas com mais de \(55\) anos que sonham em preto e branco é maior do que a proporção de pessoas com menos de \(25\) anos que sonham em preto e branco.
Passo a passo
n1 <- 306 sucessos1 <- 68 pchapeu1 <- sucessos1 / n1 cat('\npchapeu1 =', pchapeu1) ## ## pchapeu1 = 0,2222222 n2 <- 298 sucessos2 <- 13 pchapeu2 <- sucessos2 / n2 cat('\npchapeu2 =', pchapeu2) ## ## pchapeu2 = 0,04362416 dchapeu <- pchapeu1 - pchapeu2 cat('\ndchapeu =', dchapeu) ## ## dchapeu = 0,1785981 # Valor de d (diferença) na hipótese nula dzero <- 0 # Supondo a hipótese nula pbarra <- (sucessos1 + sucessos2) / (n1 + n2) erro_padrao <- sqrt( (pbarra * (1 - pbarra) / n1) + (pbarra * (1 - pbarra) / n2) ) cat('\nEP =', erro_padrao) ## ## EP = 0,02773359 # Valor p valor_p <- pnorm( dchapeu, mean = dzero, sd = erro_padrao, lower.tail = FALSE ) cat('\nValor p =', valor_p) ## ## Valor p = 0,00000000005982533Graficamente, o valor \(p\) corresponde à área em vermelho (praticamente invisível!) abaixo da seguinte distribuição normal, com média \(0\) e desvio padrão \(\approx 0{,}028\):
Como o valor \(p\) \((\approx 0{,}00000000005983)\) é muito menor do que \(\alpha\) (\(0{,}01\)), temos evidência para rejeitar \(H_0\).
Usando a função
prop.test()prop.test( matrix( c( sucessos1, sucessos2, n1 - sucessos1, n2 - sucessos2 ), ncol = 2 ), alternative = 'greater', conf.level = .99, correct = FALSE ) ## ## 2-sample test for equality of proportions without continuity correction ## ## data: matrix(c(sucessos1, sucessos2, n1 - sucessos1, n2 - sucessos2), ncol = 2) ## X-squared = 41,471, df = 1, p-value = 0,00000000005983 ## alternative hypothesis: greater ## 99 percent confidence interval: ## 0,1168363 1,0000000 ## sample estimates: ## prop 1 prop 2 ## 0,22222222 0,04362416
8 Professores \(\times\) professoras3
Professoras mulheres recebem avaliações melhores do que professores homens?
Em uma universidade, as professoras mulheres receberam as seguintes notas dos alunos:
mulheres <- c(
4.3, 4.3, 4.4, 4.0, 3.4, 4.7, 2.9, 4.0, 4.3, 3.4, 3.4, 3.3
)Na mesma universidade, os professores homens receberam as seguintes notas:
homens <- c(
4.5, 3.7, 4.2, 3.9, 3.1, 4.0, 3.8, 3.4, 4.5, 3.8, 4.3, 4.4, 4.1, 4.2, 4.0
)Nada é sabido sobre as variâncias das duas populações.
Construa um intervalo de confiança de \(95\%\) para a diferença entre as médias.
Condições
As amostras são aleatórias, representativas da população?
Supomos que sim.
Os valores das amostras são independentes entre si?
Supomos que sim.
As amostras são menos do que \(10\%\) da população?
Sim.
As amostras vêm de populações aproximadamente normais?
Podemos fazer testes de normalidade:
mulheres %>% shapiro.test() ## ## Shapiro-Wilk normality test ## ## data: . ## W = 0,91737, p-value = 0,2648homens %>% shapiro.test() ## ## Shapiro-Wilk normality test ## ## data: . ## W = 0,94474, p-value = 0,4457Os valores \(p\) acima de \(0{,}05\) são evidência de que os dados seguem uma distribuição normal.
Passo a passo
n_1 = length(mulheres) cat('n_1 =', n_1) ## n_1 = 12 xbarra_1 <- mean(mulheres) cat('\nMédia amostral 1 =', xbarra_1) ## ## Média amostral 1 = 3,866667 s_1 <- sd(mulheres) cat('\nDesvio padrão amostral 1 =', s_1) ## ## Desvio padrão amostral 1 = 0,5630006 n_2 = length(homens) cat('\nn_2 =', n_2) ## ## n_2 = 15 xbarra_2 <- mean(homens) cat('\nMédia amostral 2 =', xbarra_2) ## ## Média amostral 2 = 3,993333 s_2 <- sd(homens) cat('\nDesvio padrão amostral 2 =', s_1) ## ## Desvio padrão amostral 2 = 0,5630006 d <- xbarra_1 - xbarra_2 cat('\nDiferença =', d) ## ## Diferença = -0,1266667 erro_padrao <- sqrt(s_1^2 / n_1 + s_2^2 / n_2) cat('\nEP =', erro_padrao) ## ## EP = 0,1919363 gl <- (s_1^2 / n_1 + s_2^2 / n_2)^2 / ( (s_1^2 / n_1)^2 / (n_1 - 1) + (s_2^2 / n_2)^2 / (n_2 - 1) ) cat('\nGraus de liberdade = ', gl) ## ## Graus de liberdade = 19,06342 nivel_confianca <- .95 valor_critico <- -qt((1 - nivel_confianca)/2, df = gl) cat('\nValor crítico =', valor_critico) ## ## Valor crítico = 2,092553 margem_erro <- erro_padrao * valor_critico cat('\nME =', margem_erro) ## ## ME = 0,4016368 intervalo <- d + c(-1, 1) * margem_erro cat('\nIC = [', intervalo[1], ';', intervalo[2], ']') ## ## IC = [ -0,5283035 ; 0,2749702 ]Graficamente, o intervalo de confiança corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(19{,}063\) graus de liberdade:
Observe que esta distribuição \(t\) corresponde à variável padronizada
\[ \frac{x - (\bar x_1 - \bar x_2)}{\sqrt{s_1^2 / n_1 + s_2^2 / n_2}} \]
Trabalhamos com esta variável padronizada porque a distribuição \(t\) sempre tem média zero.
Aqui, cada valor no eixo \(x\) corresponde ao número de desvios padrão acima ou abaixo da média.
De fato, o limite inferior da área preenchida corresponde a
(intervalo[1] - (xbarra_1 - xbarra_2)) / erro_padrao ## [1] -2,092553e o limite superior da área preenchida corresponde a
(intervalo[2] - (xbarra_1 - xbarra_2)) / erro_padrao ## [1] 2,092553Usando a função
t.test()t.test( x = mulheres, y = homens ) ## ## Welch Two Sample t-test ## ## data: mulheres and homens ## t = -0,65994, df = 19,063, p-value = 0,5172 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -0,5283035 0,2749702 ## sample estimates: ## mean of x mean of y ## 3,866667 3,993333Teste a afirmação de que as médias são diferentes. Use \(\alpha = 0{,}05\).
Como o intervalo de confiança para a diferença entre as médias, calculado no item anterior, inclui o zero, já sabemos, antes de qualquer cálculo, que a hipótese de nulidade não vai ser rejeitada; isto é, as amostras são evidência de que as médias são iguais.
Passo a passo
n_1 = length(mulheres) cat('n_1 =', n_1) ## n_1 = 12 xbarra_1 <- mean(mulheres) cat('\nMédia amostral 1 =', xbarra_1) ## ## Média amostral 1 = 3,866667 s_1 <- sd(mulheres) cat('\nDesvio padrão amostral 1 =', s_1) ## ## Desvio padrão amostral 1 = 0,5630006 n_2 = length(homens) cat('\nn_2 =', n_2) ## ## n_2 = 15 xbarra_2 <- mean(homens) cat('\nMédia amostral 2 =', xbarra_2) ## ## Média amostral 2 = 3,993333 s_2 <- sd(homens) cat('\nDesvio padrão amostral 2 =', s_1) ## ## Desvio padrão amostral 2 = 0,5630006 erro_padrao <- sqrt(s_1^2 / n_1 + s_2^2 / n_2) cat('\nEP =', erro_padrao) ## ## EP = 0,1919363 d <- xbarra_1 - xbarra_2 cat('\nDiferença =', d) ## ## Diferença = -0,1266667 # Vamos padronizar esta diferença tx <- (d - 0) / erro_padrao cat('\nDiferença padronizada =', tx) ## ## Diferença padronizada = -0,6599412 gl <- (s_1^2 / n_1 + s_2^2 / n_2)^2 / ( (s_1^2 / n_1)^2 / (n_1 - 1) + (s_2^2 / n_2)^2 / (n_2 - 1) ) cat('\nGraus de liberdade =', gl) ## ## Graus de liberdade = 19,06342 valor_p <- 2 * pt( tx, df = gl ) cat('\nValor p =', valor_p) ## ## Valor p = 0,5171864Como o valor \(p\) \((\approx 0{,}517)\) é maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(19{,}063\) graus de liberdade, que corresponde à variável padronizada
\[ \frac{x - (\bar x_1 - \bar x_2)}{\sqrt{s_1^2 / n_1 + s_2^2 / n_2}} \]
Como a hipótese alternativa é \(H_A : \mu_1 - \mu_2 \neq 0\), o teste é bilateral.
Usando a função
t.test()t.test( x = mulheres, y = homens ) ## ## Welch Two Sample t-test ## ## data: mulheres and homens ## t = -0,65994, df = 19,063, p-value = 0,5172 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -0,5283035 0,2749702 ## sample estimates: ## mean of x mean of y ## 3,866667 3,993333
9 Mulheres \(\times\) maridos4
A tabela abaixo mostra as quantidades de palavras faladas em um dia pelos maridos e mulheres em \(8\) casais:
| casal | marido | mulher |
|---|---|---|
| 1 | 15.684 | 24.625 |
| 2 | 26.429 | 13.397 |
| 3 | 1.411 | 18.338 |
| 4 | 7.771 | 17.791 |
| 5 | 18.876 | 12.964 |
| 6 | 15.477 | 16.937 |
| 7 | 14.069 | 16.255 |
| 8 | 25.835 | 18.667 |
Construa um intervalo de confiança de \(95\%\) para a diferença média entre as quantidades de palavras faladas pelos dois cônjuges de um casal.
Condições
As amostras são aleatórias, representativas da população?
Supomos que sim.
Os pares são independentes entre si?
Supomos que sim. O comportamento de um casal não afeta o comportamento dos outros.
As amostras são menos do que \(10\%\) da população?
Sim.
A amostra vem de uma população aproximadamente normal?
Podemos fazer um teste de normalidade sobre as diferenças:
(df$marido - df$mulher) %>% shapiro.test() ## ## Shapiro-Wilk normality test ## ## data: . ## W = 0,97075, p-value = 0,9038O valor \(p\) acima de \(0{,}05\) é evidência de que os dados seguem uma distribuição normal.
Passo a passo
n = nrow(df) cat('\nn =', n) ## ## n = 8 dbarra <- mean(df$marido - df$mulher) cat('\nMédia amostral das diferenças =', dbarra) ## ## Média amostral das diferenças = -1677,75 s <- sd(df$marido - df$mulher) cat('\nDesvio padrão amostral das diferenças =', s) ## ## Desvio padrão amostral das diferenças = 10052,87 erro_padrao <- s / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 3554,227 gl <- n - 1 cat('\nGraus de liberdade = ', gl) ## ## Graus de liberdade = 7 nivel_confianca <- .95 valor_critico <- -qt((1 - nivel_confianca)/2, df = gl) cat('\nValor crítico =', valor_critico) ## ## Valor crítico = 2,364624 margem_erro <- erro_padrao * valor_critico cat('\nME =', margem_erro) ## ## ME = 8404,411 intervalo <- dbarra + c(-1, 1) * margem_erro cat('\nIC = [', intervalo[1], ';', intervalo[2], ']') ## ## IC = [ -10082,16 ; 6726,661 ]Graficamente, o intervalo de confiança corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(7\) graus de liberdade:
Observe que esta distribuição \(t\) corresponde à variável padronizada
\[ \frac{x - \overline{(x_1 - x_2)}}{s / \sqrt{n}} \]
Trabalhamos com esta variável padronizada porque a distribuição \(t\) sempre tem média zero.
Aqui, cada valor no eixo \(x\) corresponde ao número de desvios padrão acima ou abaixo da média.
De fato, o limite inferior da área preenchida corresponde a
(intervalo[1] - dbarra) / erro_padrao ## [1] -2,364624e o limite superior da área preenchida corresponde a
(intervalo[2] - dbarra) / erro_padrao ## [1] 2,364624Usando a função
t.test()t.test( x = df$marido, y = df$mulher, paired = TRUE ) ## ## Paired t-test ## ## data: df$marido and df$mulher ## t = -0,47204, df = 7, p-value = 0,6513 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -10082,161 6726,661 ## sample estimates: ## mean of the differences ## -1677,75Teste a afirmação de que, em casais em geral, a mulher fala mais do que o homem. Use \(\alpha = 0{,}05\).
Como o intervalo de confiança para a média das diferenças, calculado no item anterior, inclui o zero, já sabemos, antes de qualquer cálculo, que a hipótese de nulidade não vai ser rejeitada; isto é, a amostra não é evidência de que as mulheres falam mais.
Hipóteses
Chamando de \(\mu_d\) a média das diferenças \[\text{palavras}_\text{marido} - \text{palavras}_\text{mulher}\] temos
\[ H_0 : \mu_d \geq 0 \]
\[ H_A : \mu_d < 0 \]
Passo a passo
n = nrow(df) cat('\nn =', n) ## ## n = 8 dbarra <- mean(df$marido - df$mulher) cat('\nMédia amostral das diferenças =', dbarra) ## ## Média amostral das diferenças = -1677,75 s <- sd(df$marido - df$mulher) cat('\nDesvio padrão amostral das diferenças =', s) ## ## Desvio padrão amostral das diferenças = 10052,87 erro_padrao <- s / sqrt(n) cat('\nEP =', erro_padrao) ## ## EP = 3554,227 gl <- n - 1 cat('\nGraus de liberdade = ', gl) ## ## Graus de liberdade = 7 # Vamos padronizar esta diferença tx <- (dbarra - 0) / erro_padrao cat('\nDiferença padronizada =', tx) ## ## Diferença padronizada = -0,4720436 valor_p <- pt( tx, df = gl ) cat('\nValor p =', valor_p) ## ## Valor p = 0,3256283Como o valor \(p\) \((\approx 0{,}326)\) é maior do que \(\alpha\) (\(0{,}05\)), não temos evidência para rejeitar \(H_0\).
Graficamente, o valor \(p\) corresponde à área em vermelho abaixo da seguinte distribuição \(t\), com \(7\) graus de liberdade, que corresponde à variável padronizada
\[ \frac{x - \overline{(x_1 - x_2)}}{s / \sqrt{n}} \]
Como a hipótese alternativa é \(H_A : \mu_d < 0\), o teste é unilateral à esquerda.
Usando a função
t.test()t.test( x = df$marido, y = df$mulher, alternative = 'less', paired = TRUE ) ## ## Paired t-test ## ## data: df$marido and df$mulher ## t = -0,47204, df = 7, p-value = 0,3256 ## alternative hypothesis: true difference in means is less than 0 ## 95 percent confidence interval: ## -Inf 5056,012 ## sample estimates: ## mean of the differences ## -1677,75
10 Temperatura corporal com bootstrapping
A tabela abaixo mostra as temperatura corporais (em graus Celsius) de \(32\) pessoas, medidas de manhã e de tarde:
| pessoa | manhã | tarde |
|---|---|---|
| 1 | 36,7 | 36,7 |
| 2 | 37,0 | 37,1 |
| 3 | 36,3 | 36,7 |
| 4 | 36,8 | 37,1 |
| 5 | 36,8 | 37,1 |
| 6 | 36,8 | 36,4 |
| 7 | 35,9 | 37,0 |
| 8 | 36,3 | 37,0 |
| 9 | 36,3 | 37,1 |
| 10 | 36,8 | 36,7 |
| 11 | 36,3 | 36,8 |
| 12 | 36,6 | 36,7 |
| 13 | 36,9 | 36,7 |
| 14 | 36,4 | 36,1 |
| 15 | 35,7 | 36,2 |
| 16 | 36,7 | 36,8 |
| 17 | 37,1 | 37,1 |
| 18 | 37,4 | 37,4 |
| 19 | 37,1 | 36,6 |
| 20 | 36,5 | 36,9 |
| 21 | 37,2 | 36,5 |
| 22 | 36,7 | 36,8 |
| 23 | 37,1 | 36,5 |
| 24 | 36,4 | 36,2 |
| 25 | 37,2 | 36,9 |
| 26 | 36,9 | 37,0 |
| 27 | 36,2 | 36,3 |
| 28 | 37,1 | 35,9 |
| 29 | 36,3 | 36,1 |
| 30 | 37,1 | 36,9 |
| 31 | 37,2 | 36,8 |
| 32 | 37,3 | 37,0 |
O método de bootstrapping consiste em simular muitas amostras — com reposição — a partir dos elementos desta amostra original. A partir destas muitas amostras simuladas, você pode calcular intervalos de confiança como explicado abaixo:
Calcule a diferença entre a temperatura matinal e a temperatura vespertina para cada pessoa, armazenando os resultados em um vetor
diferencas, com \(32\) valores.diferencas <- df$manhã - df$tarde diferencas ## [1] 0,0 -0,1 -0,4 -0,3 -0,3 0,4 -1,1 -0,7 -0,8 0,1 -0,5 -0,1 0,2 0,3 -0,5 -0,1 ## [17] 0,0 0,0 0,5 -0,4 0,7 -0,1 0,6 0,2 0,3 -0,1 -0,1 1,2 0,2 0,2 0,4 0,3Armazene, em uma lista, \(10\) mil amostras, cada uma com \(32\) elementos sorteados com reposição a partir do vetor
diferencas(com todos os elementos tendo a mesma probabilidade de sorteio).n_amostras <- 1e4 n <- length(diferencas) amostras <- replicate( n_amostras, sample(diferencas, n, replace = TRUE), simplify = FALSE )Vamos examinar a lista gerada:
str(amostras) ## List of 10000 ## $ : num [1:32] 0 0,3 -0,3 0,4 ... ## $ : num [1:32] 0,7 0,6 0,2 -0,1 ... ## $ : num [1:32] -0,1 0 0,2 -0,1 ... ## $ : num [1:32] 0 -0,3 0,7 -0,3 ... ## $ : num [1:32] 0,6 0 0 0 ... ## $ : num [1:32] 0,6 -0,4 0 -0,4 ... ## $ : num [1:32] 0,6 1,2 0,2 -0,4 ... ## $ : num [1:32] -1,1 -0,5 0,3 0,3 ... ## $ : num [1:32] 0,3 1,2 0,2 -0,3 ... ## $ : num [1:32] -0,4 -0,1 -0,1 0 ... ## $ : num [1:32] -0,1 -0,1 0,2 0,4 ... ## $ : num [1:32] 0,3 1,2 0,4 0,2 ... ## $ : num [1:32] -0,1 0,4 0,1 -0,5 ... ## $ : num [1:32] -0,1 -0,1 0,2 0 ... ## $ : num [1:32] 0,4 -0,1 0,2 -0,1 ... ## $ : num [1:32] 0,1 0,3 -0,1 -0,8 ... ## $ : num [1:32] 0,7 -0,5 -0,7 0 ... ## $ : num [1:32] -0,8 -0,1 -0,1 1,2 ... ## $ : num [1:32] 0,7 -0,1 0,4 1,2 ... ## $ : num [1:32] -0,1 -0,4 -0,1 -0,4 ... ## $ : num [1:32] -1,1 0,2 -0,1 0,4 ... ## $ : num [1:32] -1,1 -0,3 -0,1 0,7 ... ## $ : num [1:32] 0 -0,7 0,2 -0,4 ... ## $ : num [1:32] 0,3 -0,5 -0,1 -0,4 ... ## $ : num [1:32] -0,3 0,4 0,5 0 ... ## $ : num [1:32] -0,1 -0,7 0,1 -0,1 ... ## $ : num [1:32] 0,3 -0,3 0,3 -0,1 ... ## $ : num [1:32] 0 0,3 -0,1 -0,1 ... ## $ : num [1:32] 1,2 -0,1 1,2 -0,5 ... ## $ : num [1:32] 0,4 -0,4 -0,3 -0,3 ... ## $ : num [1:32] 0,2 1,2 0,1 -0,1 ... ## $ : num [1:32] 0,4 0,2 -0,1 -0,5 ... ## $ : num [1:32] -0,1 -0,1 -0,1 0,1 ... ## $ : num [1:32] -0,1 0 -0,8 0 ... ## $ : num [1:32] 0,1 0,3 1,2 0,3 ... ## $ : num [1:32] 0 -0,7 0,6 0,2 ... ## $ : num [1:32] -0,4 -0,8 0,3 -0,4 ... ## $ : num [1:32] 0 0,1 -0,5 -0,7 ... ## $ : num [1:32] 0 -1,1 0,2 0,3 ... ## $ : num [1:32] 0,2 -1,1 0,3 -0,1 ... ## $ : num [1:32] 0,3 -0,1 0,2 -0,4 ... ## $ : num [1:32] -0,5 -0,8 0 0,6 ... ## $ : num [1:32] 0,4 0 -0,1 0 ... ## $ : num [1:32] 0 -0,1 -1,1 0,5 ... ## $ : num [1:32] -0,1 0,4 0 -0,4 ... ## $ : num [1:32] -0,1 0,7 0,6 0,1 ... ## $ : num [1:32] 0,5 0,3 -0,4 0,3 ... ## $ : num [1:32] -0,5 0,3 0,5 -0,3 ... ## $ : num [1:32] 0,2 0 -0,1 0,3 ... ## $ : num [1:32] 0,6 -0,1 0,4 -0,7 ... ## $ : num [1:32] -1,1 0,5 0,7 -0,1 ... ## $ : num [1:32] 0,2 0,2 0 -0,1 ... ## $ : num [1:32] 0,6 -0,1 -0,5 0,4 ... ## $ : num [1:32] 1,2 0,4 0,4 -0,1 ... ## $ : num [1:32] 0,1 -0,1 0 -0,5 ... ## $ : num [1:32] -0,8 -0,1 0,5 -0,3 ... ## $ : num [1:32] 0,1 0,3 -0,3 0,7 ... ## $ : num [1:32] 0,2 0 -1,1 -0,1 ... ## $ : num [1:32] -0,5 -0,1 0,5 0,4 ... ## $ : num [1:32] 1,2 -0,4 -0,3 0,2 ... ## $ : num [1:32] -0,5 0 -0,4 1,2 ... ## $ : num [1:32] 0,4 0 -0,1 0 ... ## $ : num [1:32] -0,1 0,6 0,3 -0,5 ... ## $ : num [1:32] -0,1 0,3 -1,1 -0,1 ... ## $ : num [1:32] -0,4 -0,3 -0,3 0,6 ... ## $ : num [1:32] 1,2 -0,1 0 0,4 ... ## $ : num [1:32] 0,1 -0,3 0,5 -0,8 ... ## $ : num [1:32] -0,1 -0,1 0,2 0,4 ... ## $ : num [1:32] 0,2 0,6 -1,1 0,2 ... ## $ : num [1:32] 0,7 0,6 -0,5 0,2 ... ## $ : num [1:32] 0,3 -0,1 -0,8 0,6 ... ## $ : num [1:32] 0,2 1,2 -0,4 0 ... ## $ : num [1:32] -0,3 0,7 -1,1 0,2 ... ## $ : num [1:32] -1,1 0,2 -0,1 -0,1 ... ## $ : num [1:32] -0,5 0,3 0,1 -0,1 ... ## $ : num [1:32] -0,8 1,2 -0,1 0,3 ... ## $ : num [1:32] -0,3 0 0 0,7 ... ## $ : num [1:32] -0,3 0,6 -0,1 0,2 ... ## $ : num [1:32] -0,5 0,7 -0,1 0 ... ## $ : num [1:32] -0,1 0,2 0,3 -0,5 ... ## $ : num [1:32] 0 0,5 0 -0,1 ... ## $ : num [1:32] 0 0,2 0,2 -0,1 ... ## $ : num [1:32] 1,2 0,5 0,7 0 ... ## $ : num [1:32] 0 -0,1 -0,1 0,2 ... ## $ : num [1:32] 0,5 -0,4 -0,4 -0,1 ... ## $ : num [1:32] -0,1 0,6 -0,4 -0,8 ... ## $ : num [1:32] -0,1 1,2 0,4 0,6 ... ## $ : num [1:32] 0,2 0,2 -0,1 -0,1 ... ## $ : num [1:32] -0,4 0,3 0,4 0,1 ... ## $ : num [1:32] 0,7 -0,1 0 0,3 ... ## $ : num [1:32] 0,2 -0,1 -1,1 0,4 ... ## $ : num [1:32] -0,5 0,3 -0,4 0,4 ... ## $ : num [1:32] -0,7 -0,7 0,2 0,2 ... ## $ : num [1:32] -0,7 -0,1 0 -0,7 ... ## $ : num [1:32] 0,2 0,2 -0,5 0 ... ## $ : num [1:32] 0,1 0,3 0,3 0,4 ... ## $ : num [1:32] 0,3 0,4 0,3 -0,1 ... ## $ : num [1:32] 0 0,3 -0,5 0,6 ... ## $ : num [1:32] 0,3 0,3 0 1,2 ... ## [list output truncated]Armazene, em um vetor de \(10\) mil elementos, as médias das amostras obtidas no item anterior.
As funções da família
mapaplicam uma operação (aqui, a média) a cada elemento de um vetor, retornando uma lista. Veja?purrr::mappara mais detalhes sobre programação funcional em R.medias <- amostras %>% map_dbl(~mean(.))str(medias) ## num [1:10000] 0,03125 -0,00625 0,03125 0,00625 ...Imprima os valores deste vetor que estão no quantil \(0{,}05\) e no quantil \(0{,}95\). Estes são os extremos do intervalo de confiança para a diferença média de temperatura na população, com \(95\%\) de confiança, calculados via bootstrapping.
Uma das vantagens do bootstrapping é não precisar de qualquer conhecimento sobre a população — sua distribuição, sua variância etc.
ic_bs <- medias %>% quantile(c(.05, .95)) ic_bs ## 5% 95% ## -0,134375 0,134375Use a função
t.test()do R para construir um intervalo de confiança para a diferença média de temperatura na população, com \(95\%\) de confiança — não é necessário fazer os cálculos passo a passo.ic_teste <- t.test(df$manhã, df$tarde, paired = TRUE)$conf.int ic_teste ## [1] -0,1676132 0,1676132 ## attr(,"conf.level") ## [1] 0,95Compare com o resultado do bootstrapping.
Os centros dos intervalos são
Bootstrap: \(-0{,}00000000000000077716\)
Teste: \(-0{,}00000000000000044409\)
Os extremos de um intervalo são mais ou menos próximos dos extremos do outro intervalo.
Como o bootstrapping consiste em sortear elementos da amostra original, os resultados podem ser diferentes a cada execução, mas os resultados costumam ser próximos aos dos testes paramétricos — como
prop.test.Se tomarmos uma quantidade maior de amostras — aqui, foram \(10\) mil — a precisão melhorará mais ainda. Experimente.
Baseado em https://web.archive.org/web/20101008192855/http://www.lifeslittlemysteries.com/how-is-crowd-size-estimated--1074/↩︎
Baseado em Murzyn, E. (2008). “Do we only dream in colour? a comparison of reported dream colour in younger and older adults with different experiences of black and white media.” Consciousness and Cognition, 17(4), 1228–1237. http://dx.doi.org/10.1016/j.concog.2008.09.002↩︎
Baseado em Andrew Gelman e Jennifer Hill (2009), “Replication Data for Data Analysis Using Regression Multilevel Hierarchical Models”, http://hdl.handle.net/1902.1/10285↩︎
Baseado em Mehl, M. R., Vazire, S., Nairán Ramírez-Esparza, Slatcher, R. B., & Pennebaker, J. W. (2007). “Are Women Really More Talkative Than Men?”, Science, 317(5834), 82. http://dx.doi.org/10.1126/science.1139940↩︎