Ferramentas do usuário
02_tutoriais:tutorial9:start
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior | ||
|
02_tutoriais:tutorial9:start [2022/06/21 17:31] adalardo [Cálculo do P] |
02_tutoriais:tutorial9:start [2023/09/12 10:46] (atual) |
||
|---|---|---|---|
| Linha 25: | Linha 25: | ||
| Para esses testes iremos precisar apenas de duas instrumentações poderosas que já utilizamos em outros tópicos da linguagem R: a função ''sample()'' e controle de fluxo com ciclos de iteração usando o ''for()''. | Para esses testes iremos precisar apenas de duas instrumentações poderosas que já utilizamos em outros tópicos da linguagem R: a função ''sample()'' e controle de fluxo com ciclos de iteração usando o ''for()''. | ||
| - | =====Função sample===== | + | ===== Função sample===== |
| A função ''sample()'' amostra aleatoriamente elementos de um objeto ''x''. Se não utilizarmos nenhum argumento, a função irá embaralhar os elementos do objeto, ou seja, montar um vetor de mesmo tamanho com os elementos alocados aleatoriamente. Para montar vetores de tamanhos diferentes do original precisamos indicar o tamanho do vetor resultado com o argumento ''size''. Quando colocamos o argumento ''replace = TRUE'' os elementos do vetor ''x'' são amostrados com reposição, ou seja, podem ser amostrados mais do que uma vez, sendo que no padrão ''FALSE'', cada elemento pode ser amostrado apenas um vez. Por fim, o argumento ''prob'' recebe um vetor de valores de mesmo tamanho que ''x'' e que definem a probabilidade de amostrar cada elemento do vetor original. Por exemplo, ''prob = c(0.5, 1, 1.5, 2)'', significa que o elemento ''x[4]'' tem o dobro de probabilidade de ser amostrado do que ''x[2]'' e quatro vezes mais que o ''x[1]''. | A função ''sample()'' amostra aleatoriamente elementos de um objeto ''x''. Se não utilizarmos nenhum argumento, a função irá embaralhar os elementos do objeto, ou seja, montar um vetor de mesmo tamanho com os elementos alocados aleatoriamente. Para montar vetores de tamanhos diferentes do original precisamos indicar o tamanho do vetor resultado com o argumento ''size''. Quando colocamos o argumento ''replace = TRUE'' os elementos do vetor ''x'' são amostrados com reposição, ou seja, podem ser amostrados mais do que uma vez, sendo que no padrão ''FALSE'', cada elemento pode ser amostrado apenas um vez. Por fim, o argumento ''prob'' recebe um vetor de valores de mesmo tamanho que ''x'' e que definem a probabilidade de amostrar cada elemento do vetor original. Por exemplo, ''prob = c(0.5, 1, 1.5, 2)'', significa que o elemento ''x[4]'' tem o dobro de probabilidade de ser amostrado do que ''x[2]'' e quatro vezes mais que o ''x[1]''. | ||
| Linha 42: | Linha 42: | ||
| </code> | </code> | ||
| - | ===== Dados de mandíbula de Chacal Dourado ===== | + | ===== Revisitando o teste de hipótese ===== |
| - | Agora vamos revisitar os dados de Chacal Dourado e a pergunta se há diferença no tamanho de mandíbulas entre machos e fêmeas, onde exemplificamos o teste de hipótese no [[02_tutoriais:tutorial6:start|]]. | + | Agora vamos revisitar os dados de Chacal Dourado e a pergunta se há diferença no tamanho de mandíbulas entre machos e fêmeas, onde exemplificamos o teste de hipótese no tutorial [[02_tutoriais:tutorial6:start|]]. |
| <code rsplus> | <code rsplus> | ||
| Linha 72: | Linha 72: | ||
| Se a variação encontrada é devido à variações não relacionadas ao sexo, é possível gerar essa diferença permutando os dados. Caso isso seja verdade encontraremos frequentemente diferenças iguais ou maiores que a observada. | Se a variação encontrada é devido à variações não relacionadas ao sexo, é possível gerar essa diferença permutando os dados. Caso isso seja verdade encontraremos frequentemente diferenças iguais ou maiores que a observada. | ||
| - | ===== Permutando ===== | + | No código abaixo estamos aleatorizando o vetor ''mf'' em relação ao vetor ''sexo'' e calculando a estatística de interesse novamente a partir dessa simulação, e gerando o pesudovalor em seguida. Repetimos esse procedimento algumas vezes: |
| <code rsplus> | <code rsplus> | ||
| s1.mf=sample(mf) | s1.mf=sample(mf) | ||
| Linha 89: | Linha 90: | ||
| </code> | </code> | ||
| - | ===== Criando ciclos de eventos ===== | + | ==== Criando ciclos de eventos ==== |
| + | Para repetir esse procedimento muitas vezes utilizamos um controle de fluxo com a função ''for()'', que tem a seguinte estrutura: | ||
| - | Vamos agora criar um controle de fluxo para repetir esses cálculos muitas vezes. Para isso precisamos usar a função de controle de fluxo ''for()'', que tem a seguinte estrutura: | ||
| <WRAP center round tip 60%> | <WRAP center round tip 60%> | ||
| Linha 112: | Linha 113: | ||
| - | Antes de iniciar os ciclos de iteração é desejável criar o objeto que irá armazenar os resultados de cada ciclo. Note que no caso abaixo criamos o objeto ''result'' e incluímos na sua primeira posição o valor de diferença observada entre os tamanhos médios de mandibulas de machos e fêmeas. Note também, que a variável de iteração vai assumir os valores de 2 a 1000 nos ciclos. | + | Antes de iniciar os ciclos de iteração é desejável criar o objeto que irá armazenar os resultados de cada ciclo. Note que no caso abaixo criamos o objeto ''result'' e incluímos na sua primeira posição o valor de diferença observada entre os tamanhos médios de mandibulas de machos e fêmeas. Note também que, a variável de iteração vai assumir os valores de 2 a 1000 nos ciclos. |
| Linha 132: | Linha 133: | ||
| abline(v = result[1]*-1, col = "red") | abline(v = result[1]*-1, col = "red") | ||
| </code> | </code> | ||
| - | ===== Cálculo do P ===== | + | ==== Cálculo do p-valor ==== |
| Duas perguntas distintas podem ser colocadas nesse teste de hipótese. Se há diferença entre os tamanhos ou se um tamanho é maior (menor) que outro, como já vimos no teste de hipótese. | Duas perguntas distintas podem ser colocadas nesse teste de hipótese. Se há diferença entre os tamanhos ou se um tamanho é maior (menor) que outro, como já vimos no teste de hipótese. | ||
| Linha 156: | Linha 157: | ||
| p.uni | p.uni | ||
| </code> | </code> | ||
| + | |||
| + | ==== Simula T ==== | ||
| + | |||
| + | No tutorial de [[02_tutoriais:tutorial6:start|]] também utilizamos uma função que automatiza esse teste de hipótese por simulação chamada {{ :02_tutoriais:tutorial6:simulaT.r |simulaT}}. Para relembrar, baixe a função e refaça o teste: | ||
| + | |||
| + | |||
| + | <code rsplus> | ||
| + | x11(width = 10, height = 10) | ||
| + | source("simulaT.r") | ||
| + | simulaT(macho, femea, teste = "maior", anima = TRUE) | ||
| + | </code> | ||
| + | |||
| + | As funções não precisam ser consideradas procedimento abstrato no qual não temos acesso. Uma função similar a essa foi criada durante a aula no curso de 2012, com o código que aprenderam neste tópico. Abra o arquivo da função em um editor de texto e reconheça todos os comando que estão nas linhas de código da função. Na próxima aula iremos entender como incorporar um procedimento em uma função. O primeiro passo é saber executar o procedimento em linhas de código, como fizemos no início do tutorial. | ||
| Linha 161: | Linha 175: | ||
| ===== Bootstrap ===== | ===== Bootstrap ===== | ||
| - | Vamos agora pegar o mesmo exemplo anterior e estimar o intervalo de confiança da média dos machos do chacal dourado. | + | Bootstrap é outro método de simulação computacional para calcular a imprecisão associada a uma estimativa da população estatística. O procedimento é bastante simples, amostramos com reposição o mesmo número de elementos do vetor de dados e recalculamos a estimativa de interesse. Baseado na premissa que nossa amostra é representativa da nossa população estatística, conseguimos calcular os intervalos de confiança das estimativas. |
| + | |||
| + | |||
| + | |||
| + | |||
| + | No exemplo abaixo utilizaremos os mesmos dados anteriores para exemplificar o procedimento bootstrap para calcular o intervalo de confiança da média dos machos do chacal dourado. | ||
| Primeiro vamos ver novamente esses dados e sua média: | Primeiro vamos ver novamente esses dados e sua média: | ||
| <code rsplus> | <code rsplus> | ||
| Linha 167: | Linha 187: | ||
| macho.m | macho.m | ||
| </code> | </code> | ||
| - | Agora, partindo da premissa que esses dados representam o tamanho das mandíbulas do chacal dourado, podemos fazer uma reamostragem dos nossos dados e calcular novamente a média: | + | |
| + | Agora, fazemos um aleatorização deste vetor e calcular novamente a média: | ||
| <code rsplus> | <code rsplus> | ||
| mean(sample(macho)) | mean(sample(macho)) | ||
| </code> | </code> | ||
| + | |||
| Essa média não é diferente da anterior, porque mudar a posição dos valores não afeta a estimativa da média. No entanto, se usarmos uma reamostragem com reposição (amostrar um valor e depois retorná-lo, antes de amostrar o próximo), permite que os valores já amostrados apareçam novamente na nova amostra. Vamos fazê-lo: | Essa média não é diferente da anterior, porque mudar a posição dos valores não afeta a estimativa da média. No entanto, se usarmos uma reamostragem com reposição (amostrar um valor e depois retorná-lo, antes de amostrar o próximo), permite que os valores já amostrados apareçam novamente na nova amostra. Vamos fazê-lo: | ||
| <code rsplus> | <code rsplus> | ||
| - | smacho<-sample(macho, replace=TRUE) | + | smacho <- sample(macho, replace = TRUE) |
| mean(smacho) | mean(smacho) | ||
| - | mean(sample(macho,replace=TRUE)) | + | mean(sample(macho, replace = TRUE)) |
| - | mean(sample(macho,replace=TRUE)) | + | mean(sample(macho, replace = TRUE)) |
| </code> | </code> | ||
| - | Perceba que as últimas linhas de comando produzem valores diferentes apesar de serem as mesmas. Esse processo é similar ao que usamos para fazer amostras de uma distribuição conhecida com o //rnorm()// e //rpois()//, só que agora os valores passíveis de serem amostrados são aqueles presentes nos nossos dados. | + | Perceba que as últimas linhas de comando produzem valores diferentes apesar de serem as mesmas. Esse processo é similar ao que usamos para fazer amostras de uma distribuição conhecida com o //rnorm()// e //rpois()//, só que agora os valores passíveis de serem amostrados são apenas aqueles presentes nos nossos dados. |
| - | Se repetirmos esse procedimento muitas vezes e guardarmos os resultados de cada simulação de amostras com reposição, teremos um conjunto de valores chamados pseudo-valores que representam a distribuição do nosso parâmetro e portanto podemos calcular o intervalo de confiança que desejarmos a partir dessa distribuição. | + | Se repetirmos esse procedimento muitas vezes e guardarmos os resultados de cada simulação de amostras com reposição, teremos um conjunto de pseudo-valores que representam a distribuição do nosso parâmetro e portanto, podemos calcular o intervalo de confiança que desejarmos a partir dessa distribuição. |
| - | Como repetimos uma operação muitas vezes no R? Usando novamente os ciclos produzidos pela função //for(... in ...)//, vamos fazer então 100 simulações: | + | Como repetimos uma operação muitas vezes no R? Usando novamente os ciclos produzidos pela função ''for(... in ...)'', vamos fazer então 100 simulações: |
| <code rsplus> | <code rsplus> | ||
| - | nsim=100 | + | nsim <- 100 |
| - | resulta=rep(NA,nsim) | + | resulta <- rep(NA,nsim) |
| for(i in 1:nsim) | for(i in 1:nsim) | ||
| { | { | ||
| - | resulta[i]<-mean(sample(macho, replace=TRUE)) | + | resulta[i] <- mean(sample(macho, replace = TRUE)) |
| } | } | ||
| - | ## veja os valores calculados | ||
| resulta | resulta | ||
| </code> | </code> | ||
| - | Agora só falta calcular o intervalo de confiança para o limite que interessa (95%, 99%...). Vamos calcular para um intervalo de 90%. Uma forma de faze-lo é ordenando os valores e olhado quais valores estão nos extremos com 5% de cada lado. | + | |
| + | Agora precisamos calcular o intervalo de confiança, chamado **intervalo bootstrap**, para o limite que interessa (95%, 99%...). Vamos calcular para um intervalo de 90%. Uma forma de faze-lo é ordenando os valores e olhado quais valores estão nos extremos com 5% de cada lado. | ||
| <code rsplus> | <code rsplus> | ||
| sort(resulta) | sort(resulta) | ||
| Linha 197: | Linha 221: | ||
| sort(resulta)[95] ## o valore que deixa os 5 maiores de fora | sort(resulta)[95] ## o valore que deixa os 5 maiores de fora | ||
| </code> | </code> | ||
| - | Podemos também usar a função //quantile()// definindo os quantis de interesse: | + | Podemos também usar a função ''quantile()'' definindo os quantis de interesse: |
| <code rsplus> | <code rsplus> | ||
| quantile(resulta, prob=c(0.05, 0.95)) | quantile(resulta, prob=c(0.05, 0.95)) | ||
| </code> | </code> | ||
| - | + | Definitivamente, fazer só 100 simulações, não parece adequado. Existem muitos arranjos possíveis de 10 elementos reamostrados com reposição((se não estou equivocado o número de arranjos é 10¹⁰, mas como no nosso cálculo a ordem dos elementos não importa, esse número é efetivamente menor)). Refaça o código com 1000 (mil) iterações e recalcule o intervalo. | |
| - | ===== Simula T ===== | + | |
| - | + | ||
| - | No tutorial de [[02_tutoriais:tutorial6:start|]] também utilizamos uma função que automatiza esse teste de hipótese por simulação chamada {{ :02_tutoriais:tutorial6:simulaT.r |simulaT}}. Para relembrar baixe a função e refaça o teste: | + | |
| - | + | ||
| - | + | ||
| - | <code rsplus> | + | |
| - | x11(width = 10, height = 10) | + | |
| - | source("simulaT.r") | + | |
| - | simulaT(macho, femea, teste = "maior", anima = TRUE) | + | |
| - | </code> | + | |
| - | + | ||
| - | As funções não precisam ser consideradas procedimento abstrato no qual não temos acesso. Uma função similar a essa foi criada durante a aula no curso de 2012. Abra o arquivo da função em um editor de texto e reconheça todos os comando que estão nas linhas de código da função. Na próxima aula iremos entender como incorporar um procedimento em uma função. O primeiro passo é saber executar o procedimento em linhas de código, como fizemos no início do tutorial. | + | |
| Linha 230: | Linha 242: | ||
| ===== Tesourinha e a deriva continental ===== | ===== Tesourinha e a deriva continental ===== | ||
| - | Vamos agora reproduzir a análise principal do estudo publicado na Nature em 1966 (//Geographical Distribution of the Dermaptera and the Continental Drift Hypothesis//) e descrita no primeiro capítulo do [[03_apostila:08-simulacao#fnt__1|livro do Manly]] sobre permutação. | + | |
| + | Para fechar nosso tutorial vamos reproduzir uma análise mais complexa, que foi publicada em um artigo na Nature em 1966 (//Geographical Distribution of the Dermaptera and the Continental Drift Hypothesis//) e descrita no primeiro capítulo do livro do Manly (1997 ((Manly B. F. J., 1997 Randomization, bootstrap and Monte Carlo methods in biology. 2nd Ed., Chapman and Hall, London))) sobre permutação. | ||
| A ideia era verificar se a ocorrência de tesourinhas (//Dermaptera//) estava mais correlacionada com a distribuição dos continentes atual ou antes da deriva continental. | A ideia era verificar se a ocorrência de tesourinhas (//Dermaptera//) estava mais correlacionada com a distribuição dos continentes atual ou antes da deriva continental. | ||
| A informação de interesse é a correlação da ocorrência de tesourinha entre diferentes regiões biogeográficas: Eurasia, África, Madagascar, Oriente, Austrália, Nova Zelândia, América do Sul e América do Norte. Valores positivos próximos a 1 representam composições de comunidades muito parecidas, valores próximos a -1 representam composição muito distintas. Vamos reconstruir essa matriz no objeto ''data.coef'': | A informação de interesse é a correlação da ocorrência de tesourinha entre diferentes regiões biogeográficas: Eurasia, África, Madagascar, Oriente, Austrália, Nova Zelândia, América do Sul e América do Norte. Valores positivos próximos a 1 representam composições de comunidades muito parecidas, valores próximos a -1 representam composição muito distintas. Vamos reconstruir essa matriz no objeto ''data.coef'': | ||
| Linha 257: | Linha 270: | ||
| A primeira parte da análise dos dados é calcular a correlação entre a matriz de similaridade taxonômica e de eventos de dispersão (atual e antes da deriva). | A primeira parte da análise dos dados é calcular a correlação entre a matriz de similaridade taxonômica e de eventos de dispersão (atual e antes da deriva). | ||
| - | Para isso, calculamos o coeficiente de correlação de //Pearson// entre as matrizes. Esse valor irá nos dizer se duas matrizes estão correlacionadas. A correlação pode ser positiva (até +1) se variações nos elementos de uma matriz levam a variações na mesma direção dos elementos correspondentes na outra , negativa quando em direção contrária (até -1), ou podem ser não relacionadas(0). | + | Para isso, calculamos o coeficiente de correlação de //Pearson// entre as matrizes. Esse valor irá nos dizer se duas matrizes estão correlacionadas. A correlação pode ser positiva (até +1) se variações nos elementos de uma matriz levam a variações na mesma direção dos elementos correspondentes na outra , negativa quando em direção contrária (até -1), ou podem ser não relacionadas (≅0). |
| $$ r = \frac{\sum_1^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_1^n{(x_i-\bar{x})^2}}\sqrt{\sum_1^n{(y_i-\bar{y})^2}}}$$ | $$ r = \frac{\sum_1^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_1^n{(x_i-\bar{x})^2}}\sqrt{\sum_1^n{(y_i-\bar{y})^2}}}$$ | ||
| Linha 290: | Linha 303: | ||
| cor12 ## correlação observada com a distancia atual | cor12 ## correlação observada com a distancia atual | ||
| cor13 ## correlação observada com a distancia antes da deriva | cor13 ## correlação observada com a distancia antes da deriva | ||
| - | ######################################################## | + | </code> |
| - | ### Repetir a simulação muitas vezes ################### | + | |
| - | ####################################################### | + | Para reproduzir muitas vezes o procedimento acima, vamos colocá-lo dentro de um ciclo de iteração, não sem antes criar o objeto para guardar todos os valores que queremos. |
| + | |||
| + | <code rsplus> | ||
| res.cor=data.frame(sim12=rep(NA, 5000), sim13=rep(NA,5000)) | res.cor=data.frame(sim12=rep(NA, 5000), sim13=rep(NA,5000)) | ||
| str(res.cor) | str(res.cor) | ||
| Linha 304: | Linha 320: | ||
| res.cor[s,2]<-cor(as.vector(data.sim), as.vector(dist.deriva), use="pairwise.complete.obs") | res.cor[s,2]<-cor(as.vector(data.sim), as.vector(dist.deriva), use="pairwise.complete.obs") | ||
| } | } | ||
| + | |||
| + | </code> | ||
| + | |||
| + | |||
| + | Por fim, vamos avaliar os resultados e calcular o p-valor: | ||
| + | |||
| + | <code rsplus> | ||
| str(res.cor) | str(res.cor) | ||
| par(mfrow=c(2,1)) | par(mfrow=c(2,1)) | ||
| Linha 315: | Linha 338: | ||
| p13=sum(res.cor[,2]<= res.cor[1,2])/(dim(res.cor)[1]) | p13=sum(res.cor[,2]<= res.cor[1,2])/(dim(res.cor)[1]) | ||
| p13 | p13 | ||
| - | |||
| </code> | </code> | ||
| + | |||
| + | Um fase muito importante é a interpretação dos resultados de testes como esse, que não está no escopo deste curso. De qualquer forma, consegue imaginar a conclusão do artigo para esse resultado? | ||
| ===== Para saber mais ===== | ===== Para saber mais ===== | ||
| Veja a aba da apostila deste mesmo tópico. Ali apresentamos outros conceitos. Dois livros são muito importantes e lançaram as bases das análises de Monte Carlo na ecologia: | Veja a aba da apostila deste mesmo tópico. Ali apresentamos outros conceitos. Dois livros são muito importantes e lançaram as bases das análises de Monte Carlo na ecologia: | ||
02_tutoriais/tutorial9/start.1655843470.txt.gz · Última modificação: 2022/06/21 17:31 por adalardo
Exceto onde for informado ao contrário, o conteúdo neste wiki está sob a seguinte licença: CC Attribution-Noncommercial-Share Alike 4.0 International
