Quantas amostras eu preciso coletar? E os erros tipo I e II?

Essa foi uma pergunta que eu sempre ouvi, e demorei entender.

Claro que em algumas disciplinas eu ouvi falar de número mágicos, colete 30 amostras, 10 amostras por tratamento e coisas místicas assim.

Eu também nunca entendi muito bem porque existem tantos erros, tipo I, tipo II, tipo III. É muito erro para falar se existe relação entre uma coisa e outra, a gente devia apenas estar certo ou errado. Mas com o tempo as coisas começam a fazer sentido.

Como eu so entendo (ou acho que entendo, muitas vezes estou errado) as coisas fazendo exemplos, vamos voltar as aranhas.

Vamos supor que eu queira saber se existe relação do tamanho da flor com o tamanho da aranha que fica nela. Eu poderia esperar que uma aranha grande sempre estivesse fazendo tocaia em uma flor grande enquanto aranhas pequenas estivessem sempre fazendo tocaia em flores pequenas. Vamos tentar simular esse exemplo.

Eu tenho uma medida preditora, o raio da flor, pensem em uma margarida redonda (eu sei que margaridas são varias flores, vamos abstrair isso por um momento), então para representar ela eu pego a medida o raio da circunferência que é a margarida. Depois eu vou la e pego a aranha que encontrar nela e meço o tamanho da aranha. Supondo que sempre encontre aranhas :)

Então simplificando, eu pergunto se o raio da margarida (preditor) tem uma relação linear com o tamanho (resposta) da aranhas ali de tocaia.

Então eu vou no meu campo virtual e coleto meus dados

Primeiro eu faço um gráfico pra olhar os dados:

summary(lm(resposta~preditor))

Call:
lm(formula = resposta ~ preditor)

Residuals:
     Min       1Q   Median       3Q      Max
-2.42109 -0.59796  0.05168  0.68047  2.03581

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   1.6530     0.3565   4.637 7.49e-05 ***
preditor      1.4634     0.3383   4.326 0.000174 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.052 on 28 degrees of freedom
Multiple R-squared: 0.4006,    Adjusted R-squared: 0.3792
F-statistic: 18.71 on 1 and 28 DF,  p-value: 0.0001744



E humm legal, parece que existe uma relação (eu sei, eu simulei assim).
Então eu tento criar um modelo linear:


Em vermelho temos a estatística F para o nosso modelo, com um belo p significativo e em azul temos os graus de liberdade, 1 e 28, 1 vem do fato de termos 2 parâmetros, a+bx, intercepto e inclinação, 2 parâmetros, logo 1 grau de liberdade. E denominador 28, são 30 amostras, perco um grau de liberdade por parâmetro estimado dai fica 28. Sem complicar muito aqui, o ponto é, usamos o número magico de amostras e tudo deu certo, havia um efeito, fizemos um modelo e identificamos ele.

Agora  vamos apenas simular o caso de não haver esse efeito:

  

  summary(lm(resposta.falsa~preditor))

Call:
lm(formula = resposta.falsa ~ preditor)

Residuals:
     Min       1Q   Median       3Q      Max
-2.18622 -0.61067  0.03449  0.72265  1.95235

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   2.4758     0.3493   7.089 1.03e-07 ***
preditor      0.5211     0.3314   1.572    0.127   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.031 on 28 degrees of freedom
Multiple R-squared: 0.08113,    Adjusted R-squared: 0.04831
F-statistic: 2.472 on 1 and 28 DF,  p-value: 0.1271



Tudo perfeito, notem que usamos os mesmos graus de liberdade, mas notem que agora não vemos um p significativo, o modelo não explica nada muito bem.
Ou seja representa bem o que simulamos. Parece que usando um número de amostras mágico tudo da certo.

Mas agora vamos fazer outra coisa, vamos fazer a mesma coisa que acabamos de fazer. Vamos simular apenas cinco amostras. Havia pouco tempo pra coletar, o dinheiro estava curto. Sei la, mas o que acontece com poucas amostras?

Primeiro fazendo a coleta de dados:

E ao nosso gráfico + modelo:

  summary(lm(resposta~preditor))

Call:
lm(formula = resposta ~ preditor)

Residuals:
        1         2         3         4         5
-0.080901  0.004574  0.215028  1.269947 -1.408648

Coefficients:
            Estimate Std. Error t value Pr(>|t|) 
(Intercept)   3.1261     1.0216   3.060    0.055 .
preditor     -0.6745     1.5617  -0.432    0.695 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.103 on 3 degrees of freedom
Multiple R-squared: 0.05855,    Adjusted R-squared: -0.2553
F-statistic: 0.1866 on 1 and 3 DF,  p-value: 0.6949




Primeiro olhando o grau de liberdade, vemos 1 e 3 agora, ou seja continuamos com o modelo linear de intercepto mais inclinação, mas temos 5 amostras apenas agora. Mas agora o p não é significativo. Hummmm. então simulamos um efeito, mas não conseguimos detectar ele. Ou seja nossa conclusão esta errada.
Mas e se não houve-se um efeito?
Vamos coletar os dados e olhar:

  summary(lm(resposta.falsa~preditor))

Call:
lm(formula = resposta.falsa ~ preditor)

Residuals:
      1       2       3       4       5
 0.3306  0.4478 -0.5510  0.1596 -0.3870

Coefficients:
            Estimate Std. Error t value Pr(>|t|) 
(Intercept)   1.7199     0.4749   3.622   0.0362 *
preditor      1.4490     0.7259   1.996   0.1399 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5127 on 3 degrees of freedom
Multiple R-squared: 0.5704,    Adjusted R-squared: 0.4272
F-statistic: 3.984 on 1 and 3 DF,  p-value: 0.1399


Aqui, ainda as 5 amostras, mas tenho certeza que se mostra-se esse grafico para alguem, a pessoa falaria, "Olha, parece que tem uma tendência, tem um p de 0.13, quem sabe se você coletar mais...". Mas aqui não existe um efeito, mas é fácil ser enganado e gastar tempo coletando mais para ganhar uma estrelinha (um p significativo), quem não gosta de ser premiado.

Mas chegamos ao ponto. Quanto tínhamos o número magico de amostras, nossas conclusões foram condizentes com o esperado, com aquilo que simulamos, agora com um número baixo de amostras, somos induzidos ao erro, e a dois tipos de erro. Tanto com como sem efeito, acabamos com uma tendência de concluir o contrario.

O que aconte é que:

Conclusão		Verdade
		Com Efeito	Sem Efeito
	Com Efeito	Correto	Erro Tipo I
	Sem Efeito	Erro Tipo II	Correto

Ou seja, nos tiramos uma conclusão, mas existe uma verdade que não vemos.
Quando concluimos que existe um efeito, e essa é a verdade, existe esse efeito, tudo blz, estamos corretos. Mas pode acontecer, como quando coletamos apenas 5 amostras, que existe um efeito, mas concluímos incorretamente que ele nao existe, isso é o tal do erro tipo I.

Mas a outra situação é quanto não existe um efeito, se ele realmente não existe, tudo blz, estamos corretos, mas como a gente viu aqui, podemos pensar que existe um efeito quanto ele não existe, e desperdiçar tempo atrás de mais dados inutilmente.

Mas como a gente viu, poucas amostras são ruim, muitas amostras bom, mas como se chega a esses números?
Bem existe algo que se chama "power analisys", ou a analise do poder do teste.
Melhores informações podem ser vistas aqui:
http://en.wikipedia.org/wiki/Statistical_power

Mas basicamente podemos ver qual o poder do teste que usamos. Pelo que vemos com 5 amostras parece que o teste teve um desempenho ruim, ja que nos levou ao erro, mas com 30 amostras tivemos um desempenho bom, concluímos o mesmo que simulamos.
Com o pacote pwr do R podemos ver o poder dos teste que fizemos, fornecendo o nível de significância que assumimos (0,05 normalmente), os graus de liberdade que vimos em azul e o e nível do efeito, que sempre foi 1, era o numero que multiplicava o preditor quando geramos os dados.
Nesse caso temos:

#
library(pwr)
pwr.f2.test(1,28,1,sig.level=0.05,)
pwr.f2.test(1,3,1,sig.level=0.05,)
#
Para os primeiro exemplos, notem que o u e o v, são os graus de liberdade, que pintamos de azul ali em cima, f2 é o tamanho do efeito, que usamos 1 na simulação e o sig.level é o nivel de significância, preenchendo todas essas coisas, temos qual o poder do teste.
  pwr.f2.test(1,28,1,sig.level=0.05,)

     Multiple regression power calculation

              u = 1
              v = 28
             f2 = 1
      sig.level = 0.05
          power = 0.9995525



Para o primeiros 2 exemplos, com 30 amostras, um poder de 0.999 (Máximo é 1).
Muito bom, agora quando usamos somente 5 amostras:

  pwr.f2.test(1,3,1,sig.level=0.05,)

     Multiple regression power calculation

              u = 1
              v = 3
             f2 = 1
      sig.level = 0.05
          power = 0.3434665

Apenas 0.343, bem baixo, tanto que nossas conclusões estavam saindo tudo errado. Sempre soubemos que poucas amostras eram ruim, mas agora vemos que podemos até quantificar o quão ruim é o poder dos modelos que fazemos.
Estamos usando uma regressão linear, mas para teste com proporções, anova, qualquer tipo de modelo podemos pensar qual o poder dele, a sensibilidade dele e se devemos levar a serio as conclusões que estamos chegando com os modelos.

Agora  para fechar vamos olhar uma coisa a mais, como o poder da analise muda de acordo com o numero de amostras e o tamanho do efeito que estamos observando.
Lembrando que o tamanho do efeito significa que uma flor pequena tem uma aranha miniatura, e uma flor grande tem uma aranha so um pouquinho maior, uma relação bem sutil, mas a relação pode ser também de que em uma flor pequena temos uma aranha miniatura, e uma flor grande temos uma aranha gigantesca, isso é um efeito forte, grande.

Então temos:

Primeiro vemos que independente da intensidade do efeito, quanto mais amostras melhor. Mas há um limite de bom. Chega uma hora que nem precisaríamos de mais coleta, ja teríamos um modelo com um poder "bom". Ou seja dali para frente deveríamos investir nosso tempo em analisar dados, escrever conclusões e não ficar coletando mais e mais.

Mas isso depende da intensidade do efeito, um efeito muito sutil, precisaríamos de uma quantidade muito muito grande de amostras pra tirar conclusões razoáveis, então nesse caso deveríamos é pensar bem em montar o experimento de um jeito mais eficiente.
No entanto para intensidades maiores, como a linha vermelha vemos que 30 amostras podem produzir modelos com uma sensibilidade ótima.
Claro que efeitos maiores precisam ainda de menos amostras para alcançar uma sensibilidade ótima para o teste, mas como dificilmente você sabe bem o que esperar de estimativa de intensidade, vale você ser conservador e coletar 30 amostras, se o efeito não muito fraco você vai chegar a uma boa conclusão.
E alias eu chutaria que o o número magico de 30 amostras deve ter vindo de uma analise desse tipo.

Olhar a documentação do pacote pwr vai revelar funções igual a essa pra calcular o poder de outros tipos de teste, e se continuarmos nessa linha acredito que encontraríamos o porque das sugestões de 10 amostras por nível de fator que tem em outros livros.

Mas o poder de um teste responde bem ao que ja imaginava intuitivamente, mais amostras é melhor, mas isso pode ser interessante para acabar com a velha discussão do número de amostras, e pensar se realmente vale a pena fazer algum teste dado o número de amostras que temos a mão.

Sempre o número de amostras é uma discussão recorrente mas quase nunca eu vejo a discussão acabar aqui, no poder da analise utilizada :)

  Até a próxima.




Update - Agora que me liguei que o blogger estraga tudo o código de R.
Preciso de uma solução mais eficiente que o blogger :(