Se alguém caiu aqui, eu estou no http://recologia.wordpress.com/ , não consegui fazer algumas coisas no blogger como digitar código dos scripts de R nem dar pau, no wordpress tem uma ferramenta mão na roda pra fazer isso, ai estou la :P
Mudei entre aspas que num fiquei nada aqui eheheh
quinta-feira, 8 de novembro de 2012
segunda-feira, 5 de novembro de 2012
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:
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 :(
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ãoVerdadeCom EfeitoSem EfeitoCom EfeitoCorretoErro Tipo ISem EfeitoErro Tipo IICorreto
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 :(
segunda-feira, 29 de outubro de 2012
Likelihood ou Verossimilhança.
Essa é uma palavra muito comum nos "outputs" de programas estatísticos, mas que eu demorei um pouco pra entender, se é que compreendo realmente.
Fazendo a disciplina oferecida pelo Brian Caffo no coursera (https://www.coursera.org/course/biostats) ele da um exemplo bem simples e legal.
Então vamos pensar pensar em uma moeda... na verdade não ja que moedas e bolinhas coloridas em caixinhas fazem muito pouco sentido pra mim :)
Voltando aquelas aranhas do post anterior, falamos da teoria de meta-populações, que é a forma que acho mais legal na biologia de "jogar moedas", a chance de um lugar estar ocupado ou não.
Vamos pensar que fomos olhar 4 arbustos, fomos la batemos com um pauzinho nele e em 3 arbustos caiu aranhas. Hummm
Então no primeiro estava ocupado por aranhas, o segundo idem, o terceiro não tinha aranhas e o quarto estava ocupado. Ou seja 3 arbustos com aranhas e 1 sem.
Então se pensarmos que cada arbusto tem uma chance de ter aranhas, tem sua chance de ocupação, qual a taxa de ocupação nesse exemplo?
Podemos perguntar da seguinte forma, se a taxa de ocupação for 50%, poderíamos ver esse resultado?
A resposta é sim.
Usando o R, podemos simular essa situação usando a função rbinom, que simula amostras tiradas de uma distribuição binomial.
Vamos simular 100 exemplo pra ver se algum fica igual a nossa coleta de dados.
#Vamos somar todos os casos que tivemos 3 arbustos com aranhas
#e 1 arbusto sem aranha, ou seja a amostra tem tamanhos 4, e queremos
#aquelas que em 3 encontramos aranhas
sum(rbinom(100,4,prob=0.5)==3)
E temos:
Temos 22 casos idênticos ao que coletamos.
Mas pera ai, 3 arbustos com aranhas, um sem, talvez a chance de ocupação seja maior que 50%, e se repetirmos essa simulação com 60% de chance de ocupação?
sum(rbinom(100,4,prob=0.6)==3)
E vemos 37, casos, olha ae muito mais casos ficam igual, estamos melhorando.
Mas podemos subir mais ainda, vamos tentar uma chance de ocupação de 80%
sum(rbinom(100,4,prob=0.8)==3)
Legal, agora vemos quase 50 casos, as vezes mais as vezes menos, claro que o número sempre varia, é uma simulação.
Mas a verosimilhança, ou likelihood, é exatamente isso, na verdade qualquer valor de chance de ocupação pode produzir resultados como os que observamos no campo, no entanto, a diferença é o quão comum veremos resultados idênticos aos que observamos no campo, na vida real para cada chance de ocupação que tentamos.
Ao realizar varias tentativas, vimos que algumas chances de ocupação dificilmente produziriam um resultado semelhante aquele que observamos, se tentarmos chances de ocupação ainda mais extremas, uma chance de ocupação de 20%, de 10%, ainda sim podemos observar algo como o que vimos no campo, mas é muito raro, pouco comum, então algo em torno de 80% produz mais resultados semelhantes ao que vemos.
A verosimilhança é exatamente isso, essa exploração para saber qual o valor mais provável para um parâmetro baseado num conjunto de dados, nesse caso a chance de aranhas ocuparem os arbustos, e os nossos dados são exatamente os quatro arbustos que observamos.
Então nesse caso temos, a maior verosimilhança ali perto do 80%, o pico, mas a chance de ocupação poderia assumir outros valores também, mas 20% por exemplo seria pouco provável, olhem como da um valor baixo ali no gráfico.
E 1%? Poderia também, mas a chance é quase nenhuma. Ou seja usando o likelihood da pra resumir todo o processo acima a um gráfico.
Bem, eu não vou entrar muito em formulas aqui, até porque eu sou péssimo com elas, mas da pra ter acho que ja da pra ter uma ideia melhor do que diabos é o likelihood, sendo que isso é usado pra todo tipo de parâmetro estimado.
O script para esse grafico é assim:
#
#
Copiado diretamente da aula do Brian Caffo
Mais material legal dele aqui se você pensa melhor com moedas :)
https://github.com/bcaffo/Caffo-Coursera
Fazendo a disciplina oferecida pelo Brian Caffo no coursera (https://www.coursera.org/course/biostats) ele da um exemplo bem simples e legal.
Então vamos pensar pensar em uma moeda... na verdade não ja que moedas e bolinhas coloridas em caixinhas fazem muito pouco sentido pra mim :)
Voltando aquelas aranhas do post anterior, falamos da teoria de meta-populações, que é a forma que acho mais legal na biologia de "jogar moedas", a chance de um lugar estar ocupado ou não.
Vamos pensar que fomos olhar 4 arbustos, fomos la batemos com um pauzinho nele e em 3 arbustos caiu aranhas. Hummm
Então no primeiro estava ocupado por aranhas, o segundo idem, o terceiro não tinha aranhas e o quarto estava ocupado. Ou seja 3 arbustos com aranhas e 1 sem.
Então se pensarmos que cada arbusto tem uma chance de ter aranhas, tem sua chance de ocupação, qual a taxa de ocupação nesse exemplo?
Podemos perguntar da seguinte forma, se a taxa de ocupação for 50%, poderíamos ver esse resultado?
A resposta é sim.
Usando o R, podemos simular essa situação usando a função rbinom, que simula amostras tiradas de uma distribuição binomial.
Vamos simular 100 exemplo pra ver se algum fica igual a nossa coleta de dados.
#Vamos somar todos os casos que tivemos 3 arbustos com aranhas
#e 1 arbusto sem aranha, ou seja a amostra tem tamanhos 4, e queremos
#aquelas que em 3 encontramos aranhas
sum(rbinom(100,4,prob=0.5)==3)
E temos:
Temos 22 casos idênticos ao que coletamos.
Mas pera ai, 3 arbustos com aranhas, um sem, talvez a chance de ocupação seja maior que 50%, e se repetirmos essa simulação com 60% de chance de ocupação?
sum(rbinom(100,4,prob=0.6)==3)
E vemos 37, casos, olha ae muito mais casos ficam igual, estamos melhorando.
Mas podemos subir mais ainda, vamos tentar uma chance de ocupação de 80%
sum(rbinom(100,4,prob=0.8)==3)
Legal, agora vemos quase 50 casos, as vezes mais as vezes menos, claro que o número sempre varia, é uma simulação.
Mas a verosimilhança, ou likelihood, é exatamente isso, na verdade qualquer valor de chance de ocupação pode produzir resultados como os que observamos no campo, no entanto, a diferença é o quão comum veremos resultados idênticos aos que observamos no campo, na vida real para cada chance de ocupação que tentamos.
Ao realizar varias tentativas, vimos que algumas chances de ocupação dificilmente produziriam um resultado semelhante aquele que observamos, se tentarmos chances de ocupação ainda mais extremas, uma chance de ocupação de 20%, de 10%, ainda sim podemos observar algo como o que vimos no campo, mas é muito raro, pouco comum, então algo em torno de 80% produz mais resultados semelhantes ao que vemos.
A verosimilhança é exatamente isso, essa exploração para saber qual o valor mais provável para um parâmetro baseado num conjunto de dados, nesse caso a chance de aranhas ocuparem os arbustos, e os nossos dados são exatamente os quatro arbustos que observamos.
Então nesse caso temos, a maior verosimilhança ali perto do 80%, o pico, mas a chance de ocupação poderia assumir outros valores também, mas 20% por exemplo seria pouco provável, olhem como da um valor baixo ali no gráfico.
E 1%? Poderia também, mas a chance é quase nenhuma. Ou seja usando o likelihood da pra resumir todo o processo acima a um gráfico.
Bem, eu não vou entrar muito em formulas aqui, até porque eu sou péssimo com elas, mas da pra ter acho que ja da pra ter uma ideia melhor do que diabos é o likelihood, sendo que isso é usado pra todo tipo de parâmetro estimado.
O script para esse grafico é assim:
#
thetaVals <- span="span"> seq(0, 1, length = 1000)
maxL <- span="span"> .75 ^ 3 * .25 ^1
lVals <- span="span"> thetaVals ^ 3 * (1 - thetaVals) ^1 / maxL
#postscript("coinLikelihood.eps", width=4, height=4, horizontal = FALSE)
plot(thetaVals, lVals,
type = "l",
frame = F,
xlab = expression(theta),
ylab = "Likelihood",
lwd = 2)
lines(range(thetaVals[lVals > 1 / 8]), c(1/8,1/8))
lines(range(thetaVals[lVals > 1 / 16]), c(1/16,1/16))
#
Copiado diretamente da aula do Brian Caffo
Mais material legal dele aqui se você pensa melhor com moedas :)
https://github.com/bcaffo/Caffo-Coursera
quarta-feira, 24 de outubro de 2012
Um pequeno exemplo de um modelo com Zeros Inflados
Modelos inflados por zeros podem ser muitos interessantes. Vou tentar fazer uma breve introdução com um exemplo simples.
Vamos supor a seguinte situação, eu gostaria de saber se o número de aranhas presentes numa especie de arbusto A é diferente do número de aranhas comumente encontrado num arbusto B. Utilizando o número magico de coletas que me foi passado na faculdade, eu procurei trinta arbustos da especie A e trinta arbustos da especie B e contei quantas aranhas tinham neles:
#Simulando os dados
set.seed(13)
população<-factor(rep(c("A","B"),each=30))
contagem<-rpois(60,rep(c(2,4),each=30))
população<-factor(rep(c("A","B"),each=30))
contagem<-rpois(60,rep(c(2,4),each=30))
#
Vemos que no Arbusto B, em média temos 4 aranhas, enquanto no arbusto A temos 2 aranhas.
Podemos testar essa afirmação através de um modelo de regressão de poisson.
#Modelo de Poisson
modelo1<-glm(contagem~população,family="poisson")
summary(modelo1)
exp(coef(modelo1))
summary(modelo1)
exp(coef(modelo1))
#
E obtemos como resultado:
Call:
glm(formula = contagem ~ população, family = "poisson")
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9833 -0.8255 0.0237 0.5868 2.9973
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6763 0.1302 5.195 2.05e-07 ***
populaçãoB 0.6587 0.1604 4.107 4.01e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 90.678 on 59 degrees of freedom
Residual deviance: 72.885 on 58 degrees of freedom
AIC: 235.15
Number of Fisher Scoring iterations: 5
> exp(coef(modelo1))
(Intercept) populaçãoB
1.966667 1.932203
glm(formula = contagem ~ população, family = "poisson")
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9833 -0.8255 0.0237 0.5868 2.9973
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6763 0.1302 5.195 2.05e-07 ***
populaçãoB 0.6587 0.1604 4.107 4.01e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 90.678 on 59 degrees of freedom
Residual deviance: 72.885 on 58 degrees of freedom
AIC: 235.15
Number of Fisher Scoring iterations: 5
> exp(coef(modelo1))
(Intercept) populaçãoB
1.966667 1.932203
Vejam só, o modelo nos da varias estrelinhas, nossa afirmação foi na mosca, mas mais importante que isso, temos estimativas do tamanho da população, e em média elas ficam bem próxima de realidade, o intercepto que representa a população nos arbustos A tem 1.966.. aranhas, bem próximo do 2 que usamos pra simular os dados, e no arbusto B temos 1.966 + 1.932.., que da praticamente 4, o numero que usamos na simulação. Muito legal, muitos p significativos e condecorados com variasssss estrelinhas.
Mas nem sempre o mundo é tão simples assim. Como nosso amigo Levins propôs uma vez, nem sempre todos os arbustos vão ter aranhas, eles podem não estar ocupados por aranhas.
Pra gente contar aranhas em um arbusto, tem que ter mais de uma aranha la , ou a gente vai contar zero aranhas, até ai tudo bem pois alguns arbustos podem ser um péssimo lugar pra se viver e nenhuma aranha fica la, ou seja zero aranhas.
Mas podem existir um outro tipo de zero, o zero de aranhas devido a elas nunca terem chegado la, o arbusto pode ser perfeito pra abrigar aranhas, mas nenhuma chegou la. Agora vamos colocar isso nos dados e ver o que acontece com nosso modelo.
#Simulando ocupação
ocupação<-rbinom(60,size=1,prob=0.60)
contagem.ocup<-contagem*ocupação
contagem.ocup<-contagem*ocupação
#Modelo 2 de poisson
modelo2<-glm(contagem.ocup~população,family="poisson")
summary(modelo2)
exp(coef(modelo2))
summary(modelo2)
exp(coef(modelo2))
#
Call:
glm(formula = contagem.ocup ~ população, family = "poisson")
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0656 -1.8719 -0.3133 0.6582 4.2838
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4906 0.1429 3.434 0.000594 ***
populaçãoB 0.2671 0.1898 1.407 0.159458
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 160.14 on 59 degrees of freedom
Residual deviance: 158.14 on 58 degrees of freedom
AIC: 264.01
Number of Fisher Scoring iterations: 6
> exp(coef(modelo2))
(Intercept) populaçãoB
1.633333 1.306122
glm(formula = contagem.ocup ~ população, family = "poisson")
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0656 -1.8719 -0.3133 0.6582 4.2838
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4906 0.1429 3.434 0.000594 ***
populaçãoB 0.2671 0.1898 1.407 0.159458
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 160.14 on 59 degrees of freedom
Residual deviance: 158.14 on 58 degrees of freedom
AIC: 264.01
Number of Fisher Scoring iterations: 6
> exp(coef(modelo2))
(Intercept) populaçãoB
1.633333 1.306122
É muito triste, mas o fato de nem todos os arbustos estarem ocupados por aranhas nos tira muitas muitas estrelinhas, mas pior que isso, agora concluímos que os arbustos A e B tem o mesmo numero de aranhas, isso é mal, muito mal. Olha a estimativa do intercepto 1.63, usamos 2 na simulação, mas os zeros adicionais puxaram a media do numero de aranhas nos arbustos pra baixo.
Mas o que acontece é que nesse caso não estamos levando em conta que temos 2 tipos de zeros nas coletas, existem zeros de o arbusto ter zero aranhas, mas existem zeros devido ao arbusto não estar ocupado, aranhas nunca chegaram la.
O que deveríamos fazer é tentar fazer um modelo onde esses dois tipos de situação são levados em conta.
Temos de declarar ao modelo que existem zeros devido a uma chance de ocupação que cada arbusto tem, ai considerando os arbustos ocupados que deveríamos fazer as contas para o numero de aranhas no arbusto. Essa é a solução de modelos com zero inflados tenta trazer.
#modelo zero inflado usando o pacote pscl
library(pscl)
modelo3<-zeroinfl(contagem.ocup~população|1)
summary(modelo3)
exp(coef(modelo3)[1:2])
plogis(coef(modelo3)[3])
modelo3<-zeroinfl(contagem.ocup~população|1)
summary(modelo3)
exp(coef(modelo3)[1:2])
plogis(coef(modelo3)[3])
#
Call:
zeroinfl(formula = contagem.ocup ~ população | 1)
Pearson residuals:
Min 1Q Median 3Q Max
-0.9876 -0.8934 -0.2306 0.4150 3.2704
Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.8296 0.1540 5.388 7.13e-08 ***
populaçãoB 0.5977 0.2019 2.960 0.00307 **
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4567 0.2829 -1.615 0.106
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 8
Log-likelihood: -104.9 on 3 Df
> exp(coef(modelo3)[1:2])
count_(Intercept) count_populaçãoB
2.292305 1.817889
> plogis(coef(modelo3)[3])
zero_(Intercept)
0.3877701
zeroinfl(formula = contagem.ocup ~ população | 1)
Pearson residuals:
Min 1Q Median 3Q Max
-0.9876 -0.8934 -0.2306 0.4150 3.2704
Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.8296 0.1540 5.388 7.13e-08 ***
populaçãoB 0.5977 0.2019 2.960 0.00307 **
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4567 0.2829 -1.615 0.106
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 8
Log-likelihood: -104.9 on 3 Df
> exp(coef(modelo3)[1:2])
count_(Intercept) count_populaçãoB
2.292305 1.817889
> plogis(coef(modelo3)[3])
zero_(Intercept)
0.3877701
A primeira boa noticia é que nossas estrelinhas são devolvidas, temos muitos p significativos. Mas vamos olhar as nossas estimativas, agora estimamos para o arbusto A uma média de 2.29 e para o B 2.29 + 1.81, o que é bem próximo de 4, 2 e 4 foram exatamente os valores que usamos.
Olhamos ainda a estimativa de ocupação, é de 38%, tudo bem que usamos 60%, mas é bastante razoável, e temos que levar em consideração o número de amostras, mas voltamos a concluir certo que arbustos A tem em média 2 aranhas, e arbusto B o dobro.
A parte mais interessante, é que apesar de dois processos estarem acontecendo simultaniamente, um que é o ambiente influenciando no número de aranhas que encontramos e o outro é arbustos estarem ocupados ou não, estamos modelando esses dois processos relativamente bem. Resgatando estimativas razoáveis. Se ignoramos um processo, como ignoramos o processo de ocupação no modelo 2, além de péssimas estimativas do tamanho populacional, podemos concluir erroneamente o que esta acontecendo.
Mas ao pensarmos sobre todos os processos que influenciam as nossas observações, podemos selecionar uma estrategia que leve tudo em conta, todos os processos que temos teoria para subsidiar e gerar boas estimativas.
Nesse caso, pensando que o ambiente influencia o número de aranhas que encontramos, mas encontrar aranhas também depende da ocupação, como prevê a teoria de metapopulações, uma regressão de poisson com zero inflados funciona perfeitamente. E o preço do modelo errado pode ser varias estrelinhas pior, uma conclusão errada do que está acontecendo.
Assinar:
Postagens (Atom)