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:32] adalardo [Permutando] |
02_tutoriais:tutorial9:start [2023/09/12 10:44] 127.0.0.1 edição externa |
||
---|---|---|---|
Linha 44: | Linha 44: | ||
===== Revisitando o teste de hipótese ===== | ===== 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 [[03_apostila:08-simulacao#fnt__1|livro do Manly]] 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: |