Bastão de Asclépio & Distribuição Normal

Bastão de Asclépio & Distribuição Normal

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(BSDA, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(EnvStats, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
source("eiras.friendlycolor.R")

Objetivos

Ao final desta aula (capítulo 5) o aluno deve ser capaz de:

  • discorrer sobre a lógica da testagem de hipótese nula;
  • definir hipóteses nula e alternativa;
  • definir hipótese nula para testes estatísticos bilateral e unilateral;
  • aplicar o teste estatístico de uma média populacional;
  • aplicar o teste estatístico de uma variância populacional;
  • diferenciar as probabilidades de erros do tipo I e II (alfa e beta);
  • interpretar os elementos da decisão estatística:
    • valor p;
    • nível de significância;
    • poder prospectivo;
    • tamanho de efeito populacional;
    • tamanho de amostra;
  • interpretar significância estatística e significância prática;
  • citar métodos de planejamento do estudo.

Os arquivos de dados e os códigos R podem ser utilizados explicitamente ou estão disponíveis para que os explore. Além disto, outros códigos serão criados por você e podem ser salvos com nomes que você mesmo escolherá na pasta do mesmo projeto.

Retomando a incerteza

Distribuição binomial

A distribuição binomial tem dois parâmetros, o número de tentativas (\(n\)) e a probabilidade de sucesso em cada tentativa (p).

Distribuição binomial: 1 jogada (Bernoulli)

Com apenas uma jogada temos um experimento de Bernoulli. A distribuição é composta por 50% de probabilidade de ter não ter sucesso (cara) e 50% de ter 1 sucesso (coroa).

Distribuição binomial: 2 jogadas

Repetidos experimentos de Bernoulli independentes têm distribuição binomial. Com duas jogadas há 25% de probabilidade de não haver sucesso (duas caras), 50% de probabilidade de 1 sucesso (cara e coroa ou coroa e cara) e 25% de dois sucessos (duas coroas).

Distribuição binomial: 3 jogadas

Com três jogadas há 12.5% de probabilidade de não haver sucesso (três caras), 37.5% de probabilidade de 1 sucesso (cara, cara e coroa; cara, coroa e cara; coroa, cara e cara) e 37.5% de 2 sucessos (coroa, coroa e cara; coroa, cara e coroa; cara, coroa e coroa) e 12.5% para 3 sucessos (três coroas).

A função R é dbinom(x, size, prob), indicando, respectivamente, quantas jogadas, o total de jogadas e a probabilidade de sucesso de uma jogada.

Para uma moeda balanceada (prob=0.5), a probabilidade de 0 sucesso (x=0) em 3 jogadas (size=3) é:

dbinom(x=0, size=3, prob=0.5)
[1] 0.125

1 sucesso em 3 jogadas:

dbinom(x=1, size=3, prob=0.5)
[1] 0.375

2 sucessos em 3 jogadas:

dbinom(x=2, size=3, prob=0.5)
[1] 0.375

3 sucessos em 3 jogadas:

dbinom(x=3, size=3, prob=0.5)
[1] 0.125

Mais que três jogadas não é possível: a binomial tem número máximo, finito, de jogadas. Então, para 3 jogadas, de 0 a 3 sucessos cobre toda a distribuição de probabilidades.

Distribuição binomial: 5 jogadas

Os gráficos acima foram produzidos utilizando um código R similar a este:

source("eiras.friendlycolor.R")
jogadas <- 5
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
plot(sucesso, probabilidade,
     main = paste("Binomial: ",
                  jogadas, " jogadas", sep=""),
     ylim = c(0,0.5),
     type="h", 
     col=friendlycolor(8), lwd=3)
points(sucesso,probabilidade, pch=21, 
       col=friendlycolor(8), 
       bg=friendlycolor(12))

obtendo-se:

É possível ver todos os valores em uma tabela:

jogadas <- 5
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
cat ("Sucesso\tProbabilidade\n")
for (i in 1:(jogadas+1))
{
  cat (sucesso[i],"\t",probabilidade[i],"\n")
}

obtendo-se:

Sucesso Probabilidade
0    0.03125 
1    0.15625 
2    0.3125 
3    0.3125 
4    0.15625 
5    0.03125 

… ou, mais facilmente ainda, criando um dataframe:

jogadas <- 5
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso","FR")
print(binomial)
  Sucesso      FR
1       0 0.03125
2       1 0.15625
3       2 0.31250
4       3 0.31250
5       4 0.15625
6       5 0.03125

Distribuição binomial: 15 jogadas

Este código foi modificado para usar um dataframe.

source("eiras.friendlycolor.R")
jogadas <- 15
sucesso <- 0:jogadas
probabilidade <- dbinom(sucesso,jogadas,0.5)
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso","FR")
plot(binomial$Sucesso, 
     binomial$FR,
     main = paste("Binomial: ",
                  jogadas, " jogadas", sep=""),
     ylim = c(0,0.5),
     type="h", 
     col=friendlycolor(8), lwd=3)
points(binomial$Sucesso, 
       binomial$FR,
       pch=21, 
       col=friendlycolor(8), 
       bg=friendlycolor(12))

Observe que a soma de todas as colunas é igual a 1. Então, quanto mais jogadas, maior a dispesão e menor a altura das distribuições.

(alterando a escala e os nomes dos eixos)

Com o uso de um dataframe, a função plot() joga os nomes das variáveis que nomeiam as colunas para os eixos \(x\) e \(y\). Use os parâmetros xlab e ylab:

source("eiras.friendlycolor.R")
jogadas <- 15
sucesso <- 0:jogadas
prob.sucesso <- 0.5
probabilidade <- dbinom(sucesso,jogadas,prob.sucesso)
binomial <- data.frame(sucesso,probabilidade)
names(binomial) <- c("Sucesso","FR")
plot(binomial$Sucesso, 
     binomial$FR,
     main = paste("Binomial: ",
                  jogadas, " jogadas", sep=""),
     ylim = c(0,max(binomial$FR)),
     xlab="Sucessos",
     ylab="Probabilidade",
     type="h", 
     col=friendlycolor(8), lwd=3)
points(binomial$Sucesso, 
       binomial$FR,
       pch=21, 
       col=friendlycolor(8), 
       bg=friendlycolor(12))

Distribuição binomial: 15 jogadas, moeda desbalanceada

Usando demo_15jogadas_df_p.R, que é um código R mais generalizado:

color <- 20
background <- 23
prob.sucesso <- 0.7
source("demo_15jogadas_df_p.R")

Cauda = 1

Implementamos Binomial_com_caudas.R para ilustrar a somatória das probabilidades dos valores mais extremos de uma distribuição binomial.

cauda <- 1
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0, 15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 6.103516e-05

\[\text{P}[s \le 0] + \text{P}[s \ge 15] \approx 6/100000 = 0.006\%\]

Cauda = 2

cauda <- 2
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0:1, 14:15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 0.0009765625

\[\text{P}[s \le 1] + \text{P}[s \ge 14] \approx 1/1000 = 0.1\%\]

Cauda = 3

cauda <- 3
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0:2, 13:15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 0.007385254

\[\text{P}[s \le 2] + \text{P}[s \ge 13] \approx 0.7\%\]

Cauda = 4

cauda <- 4
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0:3, 12:15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 0.03515625

\[\text{P}[s \le 3] + \text{P}[s \ge 12] \approx 3.52\%\]

Cauda = 5

cauda <- 5
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0:4, 11:15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 0.1184692

\[\text{P}[s \le 4] + \text{P}[s \ge 11] \approx 11.85\%\]

Cauda = 6

cauda <- 6
source("Binomial_com_caudas.R")

total <- sum(dbinom(c(0:5, 10:15), jogadas, p.sucesso))
cat("Total = ", total, sep="")
Total = 0.3017578

\[\text{P}[s \le 5] + \text{P}[s \ge 10] \approx 30.18\%\]

Voltando para cauda = 4

\[\text{P}[s \le 3] + \text{P}[s \ge 12] \approx 3.52\%\]

\[ \begin{cases} H_0: &\text{P}[\text{sucesso}] = 0.5 \\ H_1: &\text{P}[\text{sucesso}] \ne 0.5 \end{cases} \]

cauda <- 4
source("Binomial_com_caudas_2.R")

\(\alpha\) … probabilidade do erro do tipo I

(rejeitar \(H_0\) incorretamente)

Simulação 1

Usamos, aqui, Goodcoin.R.

source("Goodcoin.R")
Goodcoin v.20201009

Dado um valor a receber em moedas de R$1.00, metade da quantia eh oferecida em moedas com um balanceamento de referencia, e metade em moedas falsas, conhecidas por balanceamento distinto. Seu desafio eh distinguir os dois conjuntos atraves de experimentos. Numero de moedas (inteiro, default=10000): 10000

Para testar se a moeda eh verdadeira, joga-se cara ou coroa certo numero de vezes cada moeda (um experimento). Numero de lancamentos por experimento (numero inteiro, default=15): 15

Qual a proporcao maxima de moedas verdadeiras que voce aceita perder, i.e. alfa = probabilidade do erro do tipo I ou de falso-positivo). (numero entre 0 e 1). alfa (default=0.05): 0.05

As moedas verdadeiras tem balanceamento de referencia (H0). (caso queira moedas balanceadas, escolha o valor igual a 0.5) Qual a probabilidade de sortear coroa para uma moeda verdadeira? (número entre 0 e 1). P[coroa|H0] (default=0.5): 0.5

Binomial adaptada a um tratamento

Certo método educacional consegue ensinar higiene pessoal a 50% dos pacientes com autismo.

Com o novo método proposto pelo Instituto Ayres Soares, no entanto, entre 15 crianças acompanhadas, 10 (66.66%) conseguiram aprender a cuidar de sua higiene.

Se os métodos fossem iguais, somente metade das crianças (7 ou 8) deveriam aprender.

O novo método é melhor?

Verificando se o novo tratamento é melhor que o antigo

Hipótese nula e alternativa

Iniciamos presumindo que o novo tratamento é igual ao anterior:

\[ \begin{cases} H_0: &\mu_{\text{novo}} = 1 / 2 \\ H_1: &\mu_{\text{novo}} > 1 / 2 \end{cases}\\ \alpha = 0.05 \]

(este tipo de teste é chamado de unicaudal ou unilateral, pois tem direção)

alfa (\(\alpha\))

Como a distribuição binomial é discreta, nem sempre é possível conseguir a cauda com a probabilidade do erro do tipo I escolhida pelo pesquisador. Vamos encontrar a probabilidade acumulada da cauda direita abaixo de 5% para estabelecer nosso critério:

# criterio desejado
alfa <- 0.05

s <- 15+1
prob.acumulada <- 0
while (prob.acumulada < alfa)
{
  s <- s-1
  prob.acumulada <- sum(dbinom(s:15,15,0.5))
  cat("P[s>=",s,"]=",prob.acumulada,"\n",sep="")
}
s <- s+1
prob.acumulada <- sum(dbinom(s:15,15,0.5))
cat("\nCauda direita com P[s>=",s,"]=",prob.acumulada,"\n",sep="")
P[s>=15]=3.051758e-05
P[s>=14]=0.0004882813
P[s>=13]=0.003692627
P[s>=12]=0.01757812
P[s>=11]=0.05923462

Cauda direita com P[s>=12]=0.01757812

O maior valor abaixo de \(\alpha=5\%\) é associado com 12 sucessos.

Este é o valor possível nesta condição de 15 crianças:

\[\text{P}[s \ge 12] \approx 1.76\% = \alpha\]

Significa que a probabilidade de 12 ou mais crianças em 15 se beneficiarem do novo tratamento, supondo que o novo tratamento não é diferente do anterior, é menor do que \(5\%\) (ou, como é o valor possível, menor que \(1.76\%\)). Caso isto aconteça com o novo tratamento consideraremos que é uma observação improvável, pois presumimos que o novo trtamento tem o mesmo sucesso do tradicional (de 50%) e, portanto, poderemos rejeitar \(H_0\) (i.e., acreditaremos que é mais provável que o novo tratamento não seja igual ao antigo).

experimento único

Então colhemos a amostra. Apenas uma:

com este resultado:

e com esta amostra precisamos chegar à decisão estatística.

valor p

Com 10 crianças bem sucedidas, já sabemos que não será o suficiente para declarar diferença entre os tratamentos, pois não rejeitaremos \(H_0\).

Mesmo assim, podemos computar o valor p observado neste estudo: o valor p é a probabilidade de observar a melhora de 10 ou mais crianças em 15 testadas sob \(H_0\) (i.e., assumindo-se que o novo tratamento tem o mesmo efeito que o tratamento antigo):

p <- sum(binomial$FR[binomial$Sucesso>=10])
cat("p = ",p,"\n")
p =  0.1508789 

\[\text{P}[s \ge 10 | \,H_0] \approx 15.09\%\]

Em outras palavras, o valor p é uma probabilidade condicionada à hipótese nula, indicando a probabilidade de termos observado determinada diferença (neste caso, pelo menos 10 crianças) supondo que o tratamento novo é igual ao de referência (que tem sucesso em 50% das crianças e, portanto, deveria ter sucesso em 7 ou 8 crianças pelo menos metade das vezes). Então, como observar 10 ou mais crianças é provável (cerca de 15%), não acreditamos que o novo tratamento é melhor (não rejeitamos \(H_0\)) com base neste experimento.

Implementamos Binomial_H0.R para exibir as áreas de não rejeição e rejeição de \(H_0\) e o valor p:

cauda <- 4
source("Binomial_H0.R")

Em outras palavras, observar 10 crianças bem sucedidas em 15, para um tratamento que consegue sucesso em 50% das crianças, é provável (probabilidade p). Como esta probabilidade é alta (\(p>\alpha\)), não arriscaremos afirmar que estes 10 sucessos sejam devidos a um desempenho melhorado do novo tratamento:

Não se rejeita \(H_0\)

Concluindo não haver evidência de que o novo tratamento seja superior ao tratamento tradicional, tomado como referência.

Simulação 2

Supondo que o tratamento alternativo (ou uma moeda desbalanceada) tivesse 66.66% de sucesso, qual seria o comportamento da aplicação deste tratamento em 10 mil institutos, cada um deles com 15 crianças?
Goodcoin v.20201009

Dado um valor a receber em moedas de R$1.00, metade da quantia eh oferecida em moedas com um balanceamento de referencia, e metade em moedas falsas, conhecidas por balanceamento distinto. Seu desafio eh distinguir os dois conjuntos atraves de experimentos. Numero de moedas (inteiro, default=10000): 10000

Para testar se a moeda eh verdadeira, joga-se cara ou coroa certo numero de vezes cada moeda (um experimento). Numero de lancamentos por experimento (numero inteiro, default=15): 15

Qual a proporcao maxima de moedas verdadeiras que voce aceita perder, i.e. alfa = probabilidade do erro do tipo I ou de falso-positivo). (numero entre 0 e 1). alfa (default=0.05): 0.05

As moedas verdadeiras tem balanceamento de referencia (H0). (caso queira moedas balanceadas, escolha o valor igual a 0.5) Qual a probabilidade de sortear coroa para uma moeda verdadeira? (número entre 0 e 1). P[coroa|H0] (default=0.5): 0.6666

Caso a moeda fosse desbalanceada (ou se o tratamento funcionasse em cerca de \(\frac{2}{3}\) das crianças), o mais provável seria observar 10 crianças com sucesso, mas poderíamos encontrar mais ou menos: cerca de 97% das vezes observaríamos entre 6 e 13 sucessos; em outros 3% observaríamos menos que 6 ou mais que 13 sucessos.

Simulação 3

Na prática não sabemos se nosso único experimento tem sucesso em metade das vezes (é o tratamento antigo, de referência) ou em \(\frac{2}{3}\) das vezes (que é o efeito que esperamos medir, no mínimo, para declarar que o tratamento novo é melhor). Vamos colocar ambas as situações hipotéticas juntas:

source("Goodcoin.R")
Goodcoin v.20201009

Dado um valor a receber em moedas de R$1.00, metade da quantia eh oferecida em moedas com um balanceamento de referencia, e metade em moedas falsas, conhecidas por balanceamento distinto. Seu desafio eh distinguir os dois conjuntos atraves de experimentos. Numero de moedas (inteiro, default=10000): 20000

Para testar se a moeda eh verdadeira, joga-se cara ou coroa certo numero de vezes cada moeda (um experimento). Numero de lancamentos por experimento (numero inteiro, default=15): 15

Qual a proporcao maxima de moedas verdadeiras que voce aceita perder, i.e. alfa = probabilidade do erro do tipo I ou de falso-positivo). (numero entre 0 e 1). alfa (default=0.05): 0.05

As moedas verdadeiras tem balanceamento de referencia (H0). (caso queira moedas balanceadas, escolha o valor igual a 0.5) Qual a probabilidade de sortear coroa para uma moeda verdadeira? (número entre 0 e 1). P[coroa|H0] (default=0.5): 0.5

As moedas falsas tem outro balanceamento. (para simular, forneça uma probabilidade diferente de 0.5 ou deixe em branco para simular somente a moeda verdadeira) Qual a probabilidade de sortear coroa para uma moeda falsa? (número entre 0 e 1). [coroa|H1] (default=simular somente a moeda verdadeira): 0.6666

Na prática, lembramos, teremos apenas um estudo.

Como decidiremos se, por exemplo, encontrarmos 10 sucessos?

Caso o treinamento novo tenha 10 sucessos concluímos por não rejeitar a hipótese nula de que o treinamento novo é igual ao treinamento antigo, que usamos como referência. Isto é o mesmo que dizer que não encontramos evidência suficiente para dizer que o treinamento novo é melhor que o antigo.

Então, neste caso, não podemos dizer que são iguais?

Quando não rejeitamos a hipótese nula:

  • caso o treinamento novo fosse, na verdade, igual ao treinamento antigo:
    • nossa observação de 10 sucessos é um dos muitos possíveis estudos representados pela curva azul.
    • estaremos corretos ao não encontrar diferença entre os treinamentos.
    • esta probabilidade de estarmos corretos é \(1-\alpha\) (95% ou 98%).
    • em outras palavras, \(1-\alpha\) é a probabilidade de não rejeitarmos \(H_0\) corretamente.
  • caso o treinamento novo fosse, na verdade, melhor ao treinamento antigo:
    • nossa observação de 10 sucessos é um dos muitos possíveis estudos representados pela curva marrom.
    • estaremos errados, falhamos em encontrar diferença entre os treinamentos.
    • esta probabilidade de estarmos errados é \(\beta\) (80% ou 78%);
    • em outras palavras, \(\beta\) é a probabilidade de não rejeitarmos \(H_0\) incorretamente.

\(\beta\) … probabilidade do erro do tipo II

(não rejeitar \(H_0\) incorretamente)

Como decidiremos se, por exemplo, nosso estudo 12 sucessos resultasse em 12 crianças bem sucedidas com o novo treinamento?

Com 12 sucessos concluímos por rejeitar a hipótese nula de que o treinamento novo é melhor ao treinamento antigo, que usamos como referência. Isto é o mesmo que dizer que encontramos evidência suficiente para dizer que o treinamento novo é melhor que o antigo.

Podemos dizer que são diferentes (neste exemplo, que o treinamento novo é melhor)?

Quando rejeitamos a hipótese nula:

  • caso o treinamento novo fosse, na verdade, igual ao treinamento antigo:
    • nossa observação de 12 sucessos é um dos muitos possíveis estudos representados pela curva azul.
    • estaremos errados ao encontrar diferença entre os treinamentos.
    • esta probabilidade de estarmos errados é \(\alpha\) (5% ou 2%).
    • em outras palavras, \(\alpha\) é a probabilidade de rejeitarmos \(H_0\) incorretamente.
  • caso o treinamento novo fosse, na verdade, melhor ao treinamento antigo:
    • nossa observação de 12 sucessos é um dos muitos possíveis estudos representados pela curva marrom.
    • estaremos corretos em encontrar diferença entre os treinamentos.
    • esta probabilidade de estarmos corretos é \(1-\beta\) (20% ou 22%);
    • em outras palavras, \(1-\beta\) é a probabilidade de rejeitarmos \(H_0\) corretamente.

\(1-\beta\) … poder do teste

(em rejeitar \(H_0\) corretamente)
Qual é a verdade? Qual é a “realidade”?

A resposta é que nunca sabemos!

Tomada de decisão: \(\alpha\) e \(\beta\)

\(H_0\) verdadeira \(H_0\) falsa
não rejeita \(H_0\) ok \(\beta\)
rejeita \(H_0\) \(\alpha\) ok

poder (\(1-\beta\))

Probabilidade de rejeitar \(H_0\) corretamente:

O efeito
não existe
O efeito
existe
sem evidência
de efeito
não rejeitou \(H_0\)
corretamente
\(\text{confiança} = 1-\alpha\)
\(\beta\)
com evidência
de efeito
\(\alpha\) rejeitou \(H_0\)
corretamente
\(\text{poder} = 1-\beta\)

e o que acontece na prática?

Não sabemos se o efeito existe na população
não rejeitou \(H_0\),
sem evidência
de efeito, então…
… o efeito não existe e a confiança que
temos na decisão correta é \(1-\alpha\)
OU
… o efeito existe e a probabilidade
de decisão errada é \(\beta\)
(se o efeito existe -> erro do tipo II)
rejeitou \(H_0\),
com evidência
de efeito, então…
… o efeito não existe e a probabilidade
de decisão errada é \(\alpha\)
(se o efeito não existe -> erro do tipo I)
OU
… o efeito existe e a probabilidade
de decisão correta é \(1-\beta\),
i.e., o poder do teste
de declarar o efeito acertadamente.

Então, neste exemplo, como não rejeitamos \(H_0\), declarar que os dois treinamentos são iguais (i.e., aceitar \(H_0\)) é a decisão incorreta com probabilidade de …

p.sucesso <- 0.6666
jogadas <- 15
beta <- sum(dbinom(0:11, jogadas, p.sucesso))
cat("beta = ",beta,"\n")
beta =  0.7909156 

Conclusão: ESTE ESTUDO É INCONCLUSIVO.

Caso insistamos em dizer que o novo treinamento é melhor que o tradicional, temos a probabilidade \(p \approx 15\%\) de estarmos errados. Caso insistamos em dizer que os dois treinamentos têm o mesmo desempenho, temos a probabilidade \(\beta \approx 79\%\) de estarmos errados.

Não temos evidência para dizer que os dois treinamentos são diferentes, e menos ainda podemos afirmar que os dois são iguais.

O que fazer para reduzir \(\beta\)?

A literatura costuma usar o nível de significância \(\alpha=0.05=5\%\) e poder entre \(80\%\) e \(90\%\) (\(\beta\) entre \(20\%\) e \(10\%\)). Implementamos Binomial_H0H1.R para exibir a região de \(\beta\).

alfa <- 0.05
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.6666667)

ponto de corte: 11  |  12 
alfa =  5 %
beta =  79.08 %
poder =  20.92 %

Estratégia 1: aumentar o valor de \(\alpha\)

Não é boa estratégia: neste exemplo precisaríamos de valores tão altos quanto \(\alpha=0.6\)

alfa <- 0.6
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.6666667)

ponto de corte: 7  |  8 
alfa =  60 %
beta =  8.82 %
poder =  91.18 %

A probabilidade de erro do tipo I exagerada fará com que rejeitemos \(H_0\) com grande incerteza, criando a contradição de declarar a diferença entre os testes quando seus desempenhos forem numericamente iguais (com probabilidade inaceitável de errar igual a \(60\%\)).

Estratégia 2: tornar as distribuições mais estreitas

Com 78 crianças e \(\alpha=0.05\)

alfa <- 0.05
jogadas <- 78
source("Binomial_H0H1.R")

H0: Binomial(n=78,p=0.5)
H1: Binomial(n=78,p=0.6666667)

ponto de corte: 46  |  47 
alfa =  5 %
beta =  9.46 %
poder =  90.54 %

Aplicando-se o novo método aplicado a 78 crianças:

  • se 47 ou mais crianças aprenderem higiene pessoal, rejeitamos \(H_0\) e podemos afirmar que o novo treinamento é melhor que o tradicional com \(90\%\) de probabilidade de estarmos corretos (no entanto, pode ser que os treinamentos sejam iguais e estejamos errados, com \(\alpha=5\%\) de probabilidade).

  • se somente até 46 crianças aprenderem higiene pessoal, não podemos rejeitar \(H_0\). No entanto, como escolhemos poder de \(90\%\) antes do iniciarmos o estudo, podemos aceitar \(H_0\) e afirmar que os dois treinamentos são iguais porque a igualdade não deve ser decorrente de insuficiência amostral (no entanto, podemos estar enganados e, na verdade, o novo treinamento ser melhor, com \(\beta=10\%\) de probabilidade).

Notou que, agora, podemos dizer que o treinamento novo é melhor com pelo menos $${47 \over 78} = 0.6025641 \approx 60\%$$ de sucesso e antes, quando $${10 \over 15} = 0.6666667 \approx 67\%$$ das crianças tiveram sucesso, o estudo era inconclusivo?

Estratégia 3: ser capaz de detectar, somente, maiores efeitos

Caso eu não tenha como avaliar mais do que 15 crianças, o estudo só poderá distinguir o novo treinamento do tradicional caso sua performance seja ainda maior. Por tentativa e erro podemos encontrar a situação em que alcançamos pelo menos 80% de poder:

alfa <- 0.05
jogadas <- 15
pH1 <- 13/15
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.8666667)

ponto de corte: 11  |  12 
alfa =  5 %
beta =  12.92 %
poder =  87.08 %

Com apenas 15 crianças disponíveis, o estudo de um novo método poderá ser conclusivo se esperamos ter um método capaz de ensinar com sucesso \(13 / 15\) ou mais das crianças, correspondendo a poder de 87.08%. Poderá ser considerado igual ao tradicional no caso de até 11 crianças terem sucesso, ou melhor que o tradicional se 12 ou mais crianças aprenderem higiene pessoal.

Para ter uma ideia do comportamento do poder, podemos conjecturar com a variação do possível sucesso do novo treinamento:

source("Binomial_H0H1variavel.R")

Distribuição normal

A distribuição normal tem dois parâmetros, média (\(\mu\)) e desvio-padrão (\(\sigma\)).

aparência

Vamos assumir uma distribuição \(N(\mu=15,\sigma=8)\):

simetria

É uma distribuição simétrica, portanto média, moda e mediana coincidem:

Metade da área sob a curva está à esquerda e metade à direita da média.

áreas sob a curva

\(\pm 1 \text{dp}\)

cerca de 68% da área entre \(-1\) e \(+1\) desvio-padrão:

\(\pm 2 \text{dp}\)

cerca de 95% da área entre -2 e +2 desvio-padrão:

\(\pm 3 \text{dp}\)

cerca de 99.7% da área entre -3 e +3 desvio-padrão:

Por que a distribuição normal é importante em estatística?

Porque suas propriedades são conhecidas e ela aparece naturalmente sob determinadas circunstâncias. Veja, por exemplo, a aula sobre “Amostragem e Reamostragem”.

Teste de média populacional, para uma condição, bilateral

intervalos de confiança

Assuma que a distribuição da estatura do homem brasileiro em 2016 é normal com média \(\mu=173\) cm desvio-padrão \(\sigma=7\) cm.

A distribuição da estatura nesta população é:

mu <- 173
sigma <- 7
xpop <- seq(mu-4*sigma, mu+4*sigma, length.out=1000)
ypop <- dnorm(xpop, mean=mu, sd=sigma)
plot(xpop, ypop, 
     main="Distribuição normal populacional",
     xlab="Estatura (cm)",
     ylab="Densidade",
     type="l")
abline(v=mu, lty=2)

Quatro participantes de um estudo em 2020, desta mesma idade, tiveram suas estaturas medidas. Gostaríamos de saber se este grupo pode ter sido amostrado da população descrita acima quanto à média de estatura, supondo que nestes últimos quatro anos não houve alteração na estatura do homens jovens brasileiros.

Suas estaturas são 169, 174, 175 e 186 cm.

Formulamos as hipóteses nula e alternativa:

\[ \begin{cases} H_0: &\mu = 173 \\ H_1: &\mu \ne 173 \end{cases} \]

Calculamos a média amostral:

estatura <- c(169, 174, 175, 186)
n <- length(estatura)
media <- mean(estatura)
cat("média = ",media, sep="")
média = 176

A tentação é imaginar que os jovens são 3 cm mais altos atualmente, um crescimento “significativo” de 1.73% em apenas 4 anos.

No entanto, para comparar mais adequadamente a média amostral com a populacional, podemos utilizar uma métrica como o erro padrão da média (EP):

ep <- sigma / sqrt(n)
cat("ep = ",ep, sep="")
ep = 3.5

O que precisamos saber, em unidades de erro padrão, qual é a distância entre o valor populacional conhecido, \(\mu\)=173 e a média amostral media=176.

Podemos decidir pelo intervalo de confiança de 95%:

alfa <- 0.05
z <- abs(qnorm(p=alfa/2, mean=0, sd=1)) # |z| em 2.5%
IC95 = round(c(media-z*ep, media+z*ep),2)
cat("IC95 = [",IC95[1],",",IC95[2],"]", sep="")
IC95 = [169.14,182.86]

Como o valor da média populacional está dentro do intervalo de confiança, não rejeitamos a hipótese nula, portanto não temos evidência amostral de que estes 4 jovens de 19 anos não pertençam à população brasileira desta idade.

Graficamente, integrando todo o teste desenvolvido acima em demo_zbilateral1condicao.R, obtemos:

Amostra:
    média = 176
    ep = 3.5
    IC95 = [169.14,182.86]

Observe que ep é metade do valor de sigma porque a amostra tem tamanho n=4. A curva em preto é a da distribuição hipotética (analítica, paramétrica) das médias amostrais, com média igual à da amostra e desvio padrão igual a EP estimado por \(\sigma/\sqrt{n}\). Como a média populacional correspondente à \(H_0\) está dentro do IC95% estimado a partir da amostra, não rejeitamos a hipótese nula. As áreas hachuradas correspondem a 5% da área sob a curva da distribuição normal (2.5% em cada cauda). Por isso que o intervalo de confiança é de 95%.

O que aconteceria se tivéssemos uma amostra com 100 homens de 19 anos que, por acaso, tivesse a mesma média amostral de 176 cm. O código é o mesmo, exceto por utilizarmos o R para gerar a nova amostra:

set.seed(3)
estatura <- rnorm(mean=176, sd=7, n=100)

Obtendo:

Amostra:
    média = 176.0772
    ep = 0.7
    IC95 = [174.71,177.45]

É problemático representar graficamente a distribuição populacional junto à distribuição das médias amostrais, mas o fazemos aqui para que observe o comportamento da distribuição das médias amostrais com EP agora 10 vezes menor que \(\sigma\):

Portanto, a seguir representaremos somente a distribuição das médias amostrais (demo_zbilateral1condicao_100.R):

Amostra:
    média = 176.0772
    ep = 0.7
    IC95 = [174.71,177.45]

Agora o intervalo de confiança 95% não engloba a média populacional e, portanto, rejeitamos a hipótese nula de que esta amostra de 100 indivíduos veio da população geral de homens brasileiros de 19 anos.

estatística z

Tudo o que foi feito acima pode ser feito com distribuições normais padronizadas. A diferença é que usaremos a população hipotética como uma distribuição normal padronizada. Todos os valores de estatura em cm são convertidos em escores z adimensionais, e a distância em z entre a média amostral e populacional (hipotetizada) é que nos dará a decisão sobre rejeição ou não rejeição da hipótese nula.

Para a amostra com 4 indivíduos (demo_zbilateral1condicao_padrao.R) teremos:

População:
    mu = 173
    sigma = 7
Amostra:
    n = 4
    média = 176
    ep = 3.5
Distancia padronizada entre medias amostral e populacional:
    z = 0.8571429

A média amostral (m) dista 0.86 erros-padrão da média populacional (veja, no código R, a linha z <- (media-mu)/ep)). Recordando-se da anatomia da distribuição normal, 95% da área sob a normal está entre \(-1.96 \le z \le 1.96\). Então, neste caso, m está situada na região de não rejeição de \(H_0\).

Para a amostra com 100 indivíduos (demo_zbilateral1condicao_padrao100.R) teremos:

População:
    mu = 173
    sigma = 7
Amostra:
    n = 100
    média = 176
    ep = 0.7
Distancia padronizada entre medias amostral e populacional:
    z = 4.285714

Aqui, a mesma média amostral (m=176) dista 4.29 erros-padrão da média populacional. Este valor está fora do intervalo entre \(-1.96 \le z \le 1.96\) e então m está situada na região de rejeição de \(H_0\). As conclusões são as mesmas. Não rejeitamos \(H_0\) quando \(n=4\), mas rejeitamos quando \(n=100\). A diferença entre as duas situações é o valor de z, resultado da mudança de EP.

Teste z

A terceira maneira de abordar o mesmo problema é o teste z. Existe um teste z disponível em BSDA::z.test(). A novidade é o valor p.

A decisão estatística, aqui, compara o valor de p com o nível de significância \(\alpha\). As áreas hachuradas que vimos nos gráficos anteriores é de 5% (\(\alpha=5\%\)). O valor p está associado ao z calculado, correspondendo, graficamente, à área sob a curva além do valor z (área mais extrema).

O teste z foi implementado em demo_ztest.R para a amostra com 4 indivíduos:

mu <- 173
sigma <- 7
alfa <- 0.05
estatura <- c(169, 174, 175, 186)
fit <- BSDA::z.test(x=estatura, 
                    sigma.x=sigma, 
                    mu=mu,
                    alternative="two.sided", 
                    conf.level=1-alfa)
print(fit)

obtendo-se:


Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 173

Alternative Hypothesis:          True mean is not equal to 173

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 176

Data:                            estatura

Test Statistic:                  z = 0.8571429

P-value:                         0.3913659

95% Confidence Interval:         LCL = 169.1401
                                 UCL = 182.8599

Observe que o intervalo de confiança e valor de z são os mesmos calculados antes. O valor p=0.3914 é maior que \(\alpha=0.05\) e, portanto, não rejeitamos \(H_0\).

Graficamente, onde está esta área de cerca de 39.14%? Considere que o teste é bilateral e, portanto, poderíamos ter verificado a diferença entre as médias amostral e populacional como 176 - 173 = 3 ou 173 - 176 = -3. Portanto, esta área abaixo de z=-0.86 ou além de z=0.86 está assinalada em verde:

Para 100 indivíduos (demo_ztest100.R) resulta:


Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 173

Alternative Hypothesis:          True mean is not equal to 173

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 176

Data:                            estatura

Test Statistic:                  z = 4.285714

P-value:                         1.82153e-05

95% Confidence Interval:         LCL = 174.628
                                 UCL = 177.372

Aqui valor p=0.0000182153 é menor que \(\alpha=0.05\) e, portanto, rejeitamos \(H_0\).

Graficamente, as áreas abaixo de z=-4.29 ou além de z=4.29 correspondem a cerca de 0.00182% e mal são visíveis:

Novamente, as decisões com \(n=4\) e \(n=100\) foram mantidas.

Portanto, para o nível de significância \(\alpha=0.05\), as três formas de decidir são equivalentes, rejeitando-se \(H_0\) quando:

  • quando o intervalo de confiança 95% não engloba a média populacional hipotetizada;
  • pelo valor z menor que -1.96 ou maior que 1.96;
  • pelo valor \(p<\alpha\) calculado por um teste z.

Estas conclusões são referentes à significância estatística.

Em geral, estatísticas de teste tendem a aumentar (maior valor de z), correspondendo a menores valores de p quando o tamanho da amostra aumenta.

Quando \(H_0\) é rejeitada, não importa o poder a priori!

Aqui não exploramos o poder do teste. O valor de \(1-\beta\) deve ser definido durante o planejamento do experimento: é o poder prospectivo ou a priori.

Na literatura aparece o cálculo do poder retrospectivo, ou a posteriori, utilizando a estatística obtida com os dados coletados. Este cálculo do poder não faz sentido, pois será alto quando rejeitamos a hipótese nula (e, neste caso, a probabilidade de cometer erro é \(\alpha\)), e será mais baixo quando não rejeitamos a hipótese nula, que é quando deveria interessar.

O valor de \(\beta\) (recomendado ser igual ou menor que 10%), portanto, tem que ser estabelecido antes do experimento para sabermos, no caso de não rejeitarmos a hipótese nula, que a causa da não rejeição não é insuficiência amostral.

Tamanho de efeito

Além da significância estatística (existência do efeito na população), é necessário analisar, também, a significância prática (importância do efeito na população).

O tamanho de efeito pode ser visto como uma medida da magnitude do efeito que não depende do tamanho da amostra (https://en.wikipedia.org/wiki/Effect_size). Uma das medidas muito usadas é o d de Cohen. Nas palavras o próprio criador do coeficiente:

Cohen (1988)

Então, neste caso bilateral, o d de Cohen é dado por

\[d = { {|m - \mu|}\over{\sigma} }\]

ou seja, é a diferença padronizada entre as médias amostral e populacional. Obviamente, para os dois exemplos acima, o tamanho de efeito é o mesmo:

d = abs(176-173)/7
cat("d de Cohen = ",d,"\n")
d de Cohen =  0.4285714 

Adotando a seguinte tabela:

Documentação da função effectsize::interpret_cohens_d

Podemos classificar a estimativa pontual do tamanho de efeito de 0.43 como “pequeno” (\(0.20 \le d < 0.50\)). Em outras palavras, seja com 4 ou 100 indivíduos, a magnitude da diferença entre as médias amostral e populacional é pequena, do ponto de vista prático (uma diferença de 3 cm em um fenômeno que tem desvio-padrão de 7 cm é algo de pequena magnitude segundo este coeficiente). Como regra prática, se a média amostral distar mais do que \(0.5\) desvio-padrão da média populacional, a diferença começa a ter maior importância prática.

Medidas de significância prática não podem depender do tamanho da amostra.

Uma forma de obter o tamanho de efeito é eliminar a influência do tamanho da amostra do escore z dividindo-o por \(\sqrt{n}\). No exemplo das estaturas, quando \(n=4\), z=0.8571429. Então,

\(\quad d = \dfrac{|z|}{\sqrt{n}}\) = 0.4285714.

Para \(n=100\), obtemos z = 4.2857143. Então, \(d = {z/\sqrt{n}}\) = 0.4285714, que é o mesmo valor.

# n = 4
mu <- 173
sigma <- 7
alfa <- 0.05
estatura <- c(169, 174, 175, 186)
fit <- BSDA::z.test(x=estatura, 
                    sigma.x=sigma, 
                    mu=mu,
                    alternative="two.sided", 
                    conf.level=1-alfa)
print(fit)

Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 173

Alternative Hypothesis:          True mean is not equal to 173

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 176

Data:                            estatura

Test Statistic:                  z = 0.8571429

P-value:                         0.3913659

95% Confidence Interval:         LCL = 169.1401
                                 UCL = 182.8599
dC <- abs(fit$statistic)/sqrt(length(estatura))
cat("\nd de Cohen =", as.numeric(dC))

d de Cohen = 0.4285714
# n = 100
mu <- 173
sigma <- 7
alfa <- 0.05
set.seed(3)
estatura <- rnorm(mean=176, sd=sigma, n=99)
v <- 100*176-99*mean(estatura)
estatura <- c(estatura, v)
fit <- BSDA::z.test(x=estatura, 
                    sigma.x=sigma, 
                    mu=mu,
                    alternative="two.sided", 
                    conf.level=1-alfa)
print(fit)

Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 173

Alternative Hypothesis:          True mean is not equal to 173

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 176

Data:                            estatura

Test Statistic:                  z = 4.285714

P-value:                         1.82153e-05

95% Confidence Interval:         LCL = 174.628
                                 UCL = 177.372
dC <- abs(fit$statistic)/sqrt(length(estatura))
cat("\nd de Cohen =", as.numeric(dC))

d de Cohen = 0.4285714

Conceito do valor p

A decisão sobre a significância estatística utilizando o valor p é a mais habitualmente utilizada. O problema é que muitos iniciantes “robotizam” suas decisões e adotam a regra mecânica de que “se p é menor que \(5\%\), rejeita-se a hipótese nula.”

É importante ter consciência do significado do valor p:

  • constitui um dos problemas enfrentados quando conduzimos uma pesquisa o fato de não sabermos qual é o padrão existente na população de interesse.
  • o motivo de realizarmos a pesquisa é, em primeiro lugar, determinar esse padrão.
  • você precisa estar ciente de que, algumas vezes, devido ao erro amostral, obteremos padrões nas amostras que não refletem de forma acurada a população de onde as amostras foram retiradas.
  • precisamos, portanto, de um algum meio para avaliar a probabilidade de que a amostra selecionada seja um retrato fiel da população.

Os testes estatísticos nos auxiliam nesta decisão, mas isso ocorre de uma forma não de todo intuitiva.

da população para a amostra

Quando mencionamos população, as pessoas costumam ter a totalidade dos indivíduos existentes em mente. No entanto, o que define uma população depende dos critérios de inclusão e exclusão de um estudo. Por exemplo:

  • a população mundial;
  • mulheres adultas diabéticas;
  • mulheres adultas diabéticas, brasileiras, entre 30 e 50 anos de idade;
  • mulheres adultas diabéticas, brasileiras, entre 30 e 50 anos de idade, com complicações vasculares;

Também podemos ter critérios de exclusão, por exemplo, entre as mulheres selecionadas, excluímos:

  • as que têm IMC acima de 40;
  • as hipertensas;
  • as que sofreram intervenção cirúrgica nos últimos 12 meses.

Populações nem precisam ser indivíduos (humanos ou não). Por exemplo:

  • municípios brasileiros;
  • escolas estaduais de ensino médio da Cidade de São Paulo;
  • células hepáticas de camundongos.

Portanto, um conjunto de entidades de acordo com os critérios de inclusão e exclusão constituem uma população válida. Em geral, não temos acesso à população inteira; muitas vezes nem conhecemos seu perfil exato. No exemplo das mulheres diabéticas acima, no máximo teremos acesso às que sabem do diagnóstico e, na prática, àquelas cadastradas no serviço médico ao qual o pesquisador tem acesso. Precisamos, portanto, trabalhar com amostras, subconjuntos obtidos a partir da população.

A amostragem precisa ser cuidadosa, utilizando técnicas que procurem obter um subconjunto representativo da população de interesse. No entanto, por mais cuidados que tomemos, existem flutuações que podem produzir amostras que não refletem o padrão da população.

Por exemplo, observe a imagem de uma população caricata, construída artificialmente.

Esta população tem, obviamente, um padrão crescente. O que poderia acontecer com amostras obtidas desta população?

Das três amostras ilustradas aqui, as duas primeiras pegaram um padrão decrescente e um padrão neutro, que não refletem a população. A expectativa do processo de amostragem é que a probabilidade de se obter uma boa amostra (i.e., representativa da população) seja bem maior do que a de obter uma amostra que não reflita o padrão populacional.

Apresentamos, aqui, um pequeno demonstrativo que faz repetidas amostras desta população e conta quantas vezes houve sucesso ou falha em capturar um padrão crescente: experimente com demo_compadrao.R.

Na situação oposta, pode haver uma população que não tenha padrão algum, com amostras eventuais não representativas, produzindo um padrão inexistente. Experimente com demo_sempadrao.R

A replicação é uma das pedras angulares da ciência:

• Se você observa um fenômeno uma vez, então pode ter sido por acaso; se o observa duas, três ou mais vezes, pode estar começando a aprender algo sobre o fenômeno estudado
• Se o seu estudo foi o primeiro neste assunto, é sensato que você trate os resultados com certo grau de cautela
Dancey & Reidy, 2019

da amostra para a população

Como temos acesso somente às informações amostrais, o problema prático é invertido: saber se a amostra nos permite inferir sobre a população de onde, presumivelmente, a obtivemos. A certeza é impossível, então só podemos calcular probabilidades.

Um teste estatístico é o teste da hipótese nula. Assim, a primeira probabilidade a ser calculada é o valor p, que corresponde da amostra não ter padrão dado que o padrão não existe na população.

Quando existe padrão na amostra, é improvável que esta tenha vindo de uma população sem padrão, de forma que o valor p será pequeno quando condicionado à hipótese nula; a hipótese nula é rejeitada e, por exclusão, ficamos com a hipótese alternativa:

Dancey & Reidy, 2019

O contrário acontece quando a amostra não revela padrão, obtendo-se um valor p elevado; neste caso não temos evidência para rejeitar a hipótese nula:

Dancey & Reidy, 2019

efeito

O efeito pode significar a associação entre variáveis ou a diferença entre condições. A estatística tem, portanto, dois objetivos: verificar a existência do efeito e sua intensidade na população, a partir da evidência amostral, respectivamente conhecidas como significância estatística e prática.

Neste contexto, o termo mais vago, “padrão”, que utilizamos antes é substituído por efeito. Sobre os possíveis efeitos podemos definir as hipóteses.

hipóteses nula e alternativa

Todo teste estatístico enuncia, pelo menos, as hipóteses nula e alternativa. A hipótese nula ou \(H_0\) sempre declara que não existe efeito na população. A hipótese alternativa, \(H_a\) ou \(H_1\) é a nossa previsão de como condições específicas podem estar relacionadas.

Portanto, a pergunta da pesquisa costuma estar expressa em \(H_1\), mas todo o procedimento testa apenas \(H_0\). Evidência a favor de \(H_1\), consequentemente, é obtida por exclusão. Este é o motivo para que, no jargão dos estatísticos, a hipótese nula é rejeitada ou a hipótese nula não é rejeitada, e não se enuncia aceitação de nenhuma das duas.

Decisão do
pesquisador
\(H_0\): O efeito
não existe
na população
\(H_1\): O efeito
existe
na população
sem evidência
de efeito
na amostra
não rejeitou \(H_0\) corretamente
\(1-\alpha\)
(nível de confiança)
Probabilidade de erro do tipo II
\(\beta\)
com evidência
de efeito
na amostra
Probabilidade de erro do tipo I
\(\alpha\)
(nível de significância)
rejeitou \(H_0\) corretamente
\(1-\beta\)
(poder)

A regra de decisão, como vimos antes, pode ser feita pela comparação do valor p ou por intervalos de confiança. Portanto, além de formular as hipóteses, precisamos definir o valor de \(\alpha\) a priori.

Por que estabelecer \(\alpha=5\%\)?

Moore, 1997

falácias

Teste z

Retomamos o teste z bilateral utilizando as medidas de estatura (em) de uma amostra de 100 homens de 19 anos:

mu <- 173
sigma <- 7
alfa <- 0.05
set.seed(3)
estatura <- rnorm(mean=176, sd=sigma, n=99)
v <- 100*176-99*mean(estatura)
estatura <- c(estatura, v)
fit <- BSDA::z.test(x=estatura, 
                    sigma.x=sigma, 
                    mu=mu,
                    alternative="two.sided", 
                    conf.level=1-alfa)
print(fit)

Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 173

Alternative Hypothesis:          True mean is not equal to 173

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 176

Data:                            estatura

Test Statistic:                  z = 4.285714

P-value:                         1.82153e-05

95% Confidence Interval:         LCL = 174.628
                                 UCL = 177.372

Interpretando esta saída com mais detalhes:

alternative hypothesis: true mean is not equal to 173

Formamente, devemos definir as hipóteses nula e alternativa antes de executar um teste. No entanto, esta saída relata, de acordo com a forma como chamamos a função R, a hipótese alternativa “média verdadeira não é igual a 173”. Por contraste, a hipótese nula precisa ser “média verdadeira igual a 173”. O teste é bilateral Podem ser redigidas assim:

\[ \begin{cases} H_0: &\mu = 173 \\ H_1: &\mu \ne 173 \end{cases} \\ \alpha=0.05 \]

Observe que o termo “verdadeira” neste contexto significa “populacional”. As hipóteses devem ser escritas com letras gregas de acordo com a convenção que as reserva para a população. O teste é feito com a amostra, para inferir seu resultado populacional.

O nível de significância não é informado para esta função. Deve ser relatado e estar em mente do pesquisador para a tomada de decisão. O teste, como vimos, fornece os valores de intervalo de confiança 95% e de z, além do valor p associado, correspondendo à área sob a curva igual ou mais extrema que o valor z observado.

O valor p é a probabilidade de observarmos determinado afastamento (z) sob a hipótese nula \(H_0\). Isto é dizer que, se não existisse diferença entre a média populacional e a média observada na amostra, o valor z encontrado tem probabilidade p de acontecer. Caso o valor de p seja grande, maior do que o nível de significância (\(p \ge \alpha\)), este valor é provável e não há evidência para rejeitar \(H_0\). Caso contrário, se \(p < \alpha\), é improvável que a diferença seja observada, rejeitamos \(H_0\) e temos evidência que apoia a existência desta diferença.

teste z unilateral para uma condição

Quando a direção do relacionamento ou da diferença (efeito) é especificada, então o teste é unilateral/unicaudal; caso contrário, é bilateral/bicaudal. Em geral (mas nem sempre), se você tiver obtido um valor p para um teste bilateral e quiser saber o valor mínimo correspondente para o teste unilateral, então: \(p_{\text{unilateral}} \ge \dfrac{p_{\text{bilateral}}}{2}\).

Exemplo:

Um pesquisador deseja investigar se o etilismo durante a prenhez diminui o peso dos ratos recém-nascidos (PRRN).

PRRN tem distribuição normal com média e o desvio-padrão populacionais sem etilismo na prenhez iguais, respectivamente, a 20 g e 4 g. Essas estimativas foram obtidas de uma amostra de 20 PRNs independentes.

50 fêmeas foram emprenhadas num experimento no qual elas receberam doses diárias de álcool. Um rato recém-nascido de cada fêmea foi selecionado aleatoriamente para compor a amostra de 50 PRNs independentes.

A média do PRRN das fêmeas que receberam doses diárias de álcool no experimento é 18 g.

Teste

  • Média populacional do PRRN hipotetizada: \(\mu = 20\)

Suposições

  • PRRN tem distribuição normal (desnecessária, devido ao teorema central do limite, TCL, pois \(n > 30\))
  • Desvio-padrão \(\sigma=4\) conhecido
  • \(n = 50\) observações independentes
  • Teste unilateral
  • Nível de confiança de 95% (ou nível de significância adotado de 5%)

Hipóteses

\[ \begin{cases} H_0: &\mu = 20 \\ H_1: &\mu < 20\;\text{(baixo PRRN)} \end{cases} \]

Estatísticas

Implementado em demo_zunilateral_PRRN.R:

Populacao:
    mu = 20
    sigma = 4
Amostra:
    n = 50
    média = 18.004
    ep = 0.5656854
Teste z:
Results of Hypothesis Test
--------------------------

Null Hypothesis:                 mean = 20

Alternative Hypothesis:          True mean is less than 20

Test Name:                       One-sample z-Test

Estimated Parameter(s):          mean of x = 18.004

Data:                            dados

Test Statistic:                  z = -3.528463

P-value:                         0.0002089903

95% Confidence Interval:         LCL =       NA
                                 UCL = 18.93447

Significancia pratica:
    d de Cohen = -0.499 (pequeno)

Note que neste teste unilateral, o limite inferior é indefinido (só interessa o valor máximo).

Decisão

  • Critério do valor crítico da estatística de teste: Como z=-3.53 < -1.64, rejeitar \(H_0\);
  • Critério do IC95: Como IC95 [NA,18.93] não contém 20 (o valor máximo não o alcança), rejeitar \(H_0\);
  • Critério do valor p: Como o valor p unilateral (less) = 2.1e-04 é menor que 5%, rejeitar \(H_0\).

“O valor-p é a probabilidade de que a estatística de teste seja igual ou mais extrema que o valor observado na direção prevista pela hipótese alternativa (\(H_1\)), presumindo que a hipótese nula (\(H_0\)) é verdadeira.”


Agresti & Finlay, 2012

O teste é unilateral (só investigamos se há redução da peso dos ratos recém-nascidos) e a direção é explícita em \(H_1\), com o símbolo \(<\). É comum encontrar a anotação da hipótese nula com o símbolo complementar (\(\ge\) neste exemplo):

\[ \begin{cases} H_0: &\mu \ge 20 \\ H_1: &\mu < 20 \; \text{(baixo PRRN)} \end{cases} \]

Optamos por usar o símbolo de igualdade (\(=\)) porque mais adequadamente espelha o que se espera de \(H_0\), a ausência de efeito. Matematicamente, também, é equivalente (Gatás, 1978, p. 220-3)

Teste de desvio-padrão populacional

No exemplo do estudo sobre o baixo peso ao nascer de ratos submetidos a etilismo, o desvio-padrão populacional foi informado (\(\sigma=4\)). Podemos testar se, de acordo com a amostra, o desvio-padrão populacional deve ser este.

Teste

  • Desvio-padrão populacional do PRRN hipotetizada: \(\sigma = 4\)

Suposições

  • PRRN tem distribuição normal (desnecessária pelo TCL)
  • \(n = 50\) observações independentes
  • Teste bilateral
  • Nível de confiança de 95% (ou nível de significância adotado de 5%)

Hipóteses

\[ \begin{cases} H_0: &\sigma = 4 \\ H_1: &\sigma \ne 4 \end{cases} \]

Estatísticas

Implementado em demo_qui2dp_PRRN.R, utilizando o teste de variância EnvStats::varTest() adaptado:

Teste de variancia:

Results of Hypothesis Test
--------------------------

Null Hypothesis:                 variance = 16

Alternative Hypothesis:          True variance is not equal to 16

Test Name:                       Chi-Squared Test on Variance

Estimated Parameter(s):          variance = 26.13468

Data:                            dados

Test Statistic:                  Chi-Squared = 80.03745

Test Statistic Parameter:        df = 49

P-value:                         0.006746438

95% Confidence Interval:         LCL = 18.23633
                                 UCL = 40.58319

Conversão em desvio-padrão:
    Desvio-padrao amostral = 5.112209
    IC95(dp) = [4.270402,6.370494]

Decisão

Os três critérios podem ser usados:

  • Critério do valor crítico da estatística de teste utiliza a distribuição qui-quadrado (\(\chi^2\)), que será visto em capítulos adiante e vamos dispensá-lo neste momento (os outros critérios bastam, são equivalentes);
  • Critério do IC95: como IC95(dp) [4.27,6.37] não contém 4, rejeitar \(H_0\);
  • Critério do valor p: Como o valor p bilateral = 6.75e-03 é menor que 5%, rejeitar \(H_0\).

Planejamento do estudo em G*Power

Para o planejamento existe um programa chamado G*Power.

Habitualmente as pessoas pensam apenas em calcular o tamanho da amostra após escolherem os valores de \(\alpha\) e \(\beta\). Comentamos acima, quando apresentávamos a distribuição binomial, situações em que o pesquisador tem determinado número de pacientes possíveis, e nada adianta um estatístico impor outro tamanho de amostra.

Quatro elementos entram no cálculo amostral: nível de significância (\(\alpha\)), poder do teste (1-\(\beta\)), significância prática (tamanho de efeito a ser detectado) e tamanho da amostra (\(n\)). Dados quaisquer três, o quarto é determinado. Estão errados os pesquisadores que pensam que o cálculo da amostra é o cálculo do tamanho da amostra. A tabela seguinte mostra outras possibilidades, todas oferecidas pelo G*Power:

Coelho, 2008

Vamos exercitar o G*Power nas mesmas situações que mostramos com o exemplo das crianças autistas aprendendo higiene pessoal.

\(n=15\), detectando efeito de 2/3 em comparação com 1/2

Computamos, nesta condição, que \(\beta\) era alto e, consequentemente, seu complemento - o poder do teste - era baixo.

pH0 <- 0.5
pH1 <- 2/3
alfa <- 0.05
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.6666667)

ponto de corte: 11  |  12 
alfa =  5 %
beta =  79.08 %
poder =  20.92 %

Imaginando que estávamos na fase de planejamento, para verificar que o poder seria insuficiente no G*Power:

  • indicar o tipo de teste: “Proportion: Difference from constant (binomial test, one sample case)”.
  • escolher a estratégia: “Compute achieved power - given \(\alpha\), sample size, and effect size”,
  • informar que é unilateral: “Tails(s) = One”,
  • tamanho do efeito: Effect size g = 0.16666 (pois buscamos 2/3 - 1/2),
  • nível de significância: “\(\alpha~\text{err prob}\) = 0.05”
  • \(n\): “Total sample size = 15”

A clicar o botão [Calculate], o G*Power mostra no lado direito da tela o ponto de corte, o poder e o \(\alpha\) possível (\(\approx 1.76\%\)), com os mesmos valores que calculamos antes.

aumentado \(\alpha\)

Buscando aumentar o poder, tentamos algumas estratégias. A primeira era aumentar o valor de \(\alpha\).

pH0 <- 0.5
pH1 <- 2/3
alfa <- 0.6
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.6666667)

ponto de corte: 7  |  8 
alfa =  60 %
beta =  8.82 %
poder =  91.18 %

O poder chega aos 90%, mas o valor de \(\alpha\) é tão alto que o novo treinamento, atingindo o mesmo desempenho médio do treinamento de referência (8 crianças com sucesso em 15), rejeitaríamos \(H_0\) com o valor p de 50% - não faz sentido!

No G*Power, existe a possibilidade de perguntar qual o \(\alpha\) a ser utilizado, dado o poder, tamanho de efeito e tamanho da amostra (estratégia que não costuma ser vista na literatura):

Obtivemos a mesma resposta: caso o número de crianças disponível fosse 15, sem possibilidade de aumentar a amostra (o que pode acontecer, por exemplo, com uma doença pouco frequente), saberíamos que o teste teria baixo poder (como antes) ou teria que usar um valor de \(\alpha\) inaceitável; nestas condições, nem deveríamos perder tempo tentando realizar este estudo.

aumentando o tamanho da amostra

Na sequência buscamos aumentar o tamanho da amostra. Por tentativa-e-erro vimos que com \(n=78\) conseguíamos, mantendo \(\alpha=0.05\), o poder mínimo de 90%.

pH0 <- 0.5
pH1 <- 2/3
alfa <- 0.05
jogadas <- 78
source("Binomial_H0H1.R")

H0: Binomial(n=78,p=0.5)
H1: Binomial(n=78,p=0.6666667)

ponto de corte: 46  |  47 
alfa =  5 %
beta =  9.46 %
poder =  90.54 %

No G*Power é mais fácil. Basta solicitar “A priori: Compute required sample size - given \(\alpha\), power, and effect size”. Esta é a abordagem que a maioria dos pesquisadores adota:

Observe o tamanho da amostra calculada no lado direito da tela: 78 participantes é a quantidade adequada para atender \(\alpha=0.05\), \(\text{poder}=1-\beta=0.90\), tamanho de efeito (\(g=0.166666\approx {2 \over 3}-{1 \over 2}\)).

incapacidade de aumentar o tamanho da amostra

Nossa terceira opção foi conformar-se em saber qual é o tamanho de efeito que seríamos capazes de detectar, dado que não podemos dispor de mais que 15 crianças. Por tentativa-e-erro encontramos que éramos capazes de distinguir treinamentos que dessem 88% de sucesso em comparação com o treinamento de referência, bem sucedido em 50% das tentativas:

pH0 <- 0.5
pH1 <- 0.88
alfa <- 0.05
jogadas <- 15
source("Binomial_H0H1.R")

H0: Binomial(n=15,p=0.5)
H1: Binomial(n=15,p=0.88)

ponto de corte: 11  |  12 
alfa =  5 %
beta =  9.59 %
poder =  90.41 %

No G*Power este valor é facilmente encontrado com a escolha de “Sensitivity: Compute required effect size - given \(\alpha\), power, and sample size”.

Verifique, na saída à direita, que o G*Power informa “Effect size g = 0.3782313”, aproximadamente o mesmo que escolhemos acima (\(0.88-0.50=0.38\)). Portanto, em um estudo no qual não temos como obter mais do que 15 pacientes, já saberíamos que nada encontraríamos se o novo treinamento não tivesse pelo menos 88% de sucesso médio. Caso nossos indícios mostrassem que o efeito não seria tão grande ou se diferenças menores não nos interessassem, poderíamos evitar levar este estudo adiante.

Referências

Agresti A & Finlay B (2012) Métodos estatísticos para as Ciências Sociais. Porto Alegre: Penso.

Coelho JP et al. (2008) Inferência Estatística: com utilização do SPSS e GPower. Lisboa: Sílabo.

Cohen J (1988) Statistical power analysis for the behavioral sciences. 2nd ed. NY: IEA.

Dancey CP & Reidy J (2019) Estatística sem Matemática para Psicologia. 7ª ed. Porto Alegre: Penso.

Gatás RR (1978) Elementos de Probabilidade e Inferência. SP: Atlas.

Moore DS (1997) Statistics: Concepts and Controversies, 4th ed. NY: Freeman.

van Belle G (2008) Statistical rules of thumb. 2nd ed. NJ: Wiley