Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- 02_tutoriais:tutorial9:start [2020/08/12 06:04]
127.0.0.1 edição externa
+++ 02_tutoriais:tutorial9:start [2023/09/12 10:46] (atual)
@@ Linha 1: / Linha 1: @@
 <WRAP tabs>
-  * [[bie5782:02_tutoriais:tutorial9:start|Tutorial]]
+  * [[02_tutoriais:tutorial9:start|Tutorial]]
-  * [[bie5782:01_curso_atual:exerPermuta| Exercícios]]
+  * [[01_curso_atual:exerPermuta| Exercícios]]
-  * [[bie5782:03_apostila:08-simulacao| Apostila]]
+  * [[03_apostila:08-simulacao| Apostila]]
 </WRAP>
 ====== 8. Reamostragem e Simulação ======
-=====Função sample=====
-Criar um vetor de LETTERS com letras de "A" a "J" e aplique a ele a função ''sample''.
-Se não utilizar nenhum argumento ela apenas embaralha os atributos do objeto. Com o argumento "replace=T", ela reamostra cada elemento com reposição, e
-o argumento "prob" é a reamostragem de cada elemento com probabilidades diferentes.
-<code>
+<WRAP center round box 70%>
-vetor=rep(LETTERS[1:10])
+{{  :02_tutoriais:tutorial9:dif.gif?|  }}
+</WRAP>
+Os métodos de Monte Carlo são procedimentos de simulação computacional para soluções de problemas complexos. Eles são utilizados em principalmente três campos da matemática: otimização, integração numérica e aleatorização de amostras de funções probabilísticas. Aqui vamos focar em sua base mais simples: a aleatorização ou permutação para gerar distribuições probabilísticas do cenário nulo e calcular a probabilidade do resultado obtido ter sido gerado pelo acaso e/ou calcular o intervalos de confiança de alguma estatística de interesse.
+Especificamente para o testes de hipótese por aleatorização ou randomização dos dados, temos ao menos as seguintes etapas:
+<WRAP center round box 60%>
+  - Definir a estatística de interesse (EI)
+  - Calcular a estatística de interesse a partir dos dados
+  - Estabelecer o cenário nulo
+  - Simular cenário nulo
+  - Calcular a EI no cenário nulo (pseudovalor)
+  - Produzir a distribuição dos pseudovalores
+  - Posicionar a EI observada na distribuição dos pseudovalores
+  - Calcular o p-valor
+</WRAP>
+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=====
+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]''.
+Vamos criar um vetor a partir do objeto ''LETTERS'', com letras de "A" a "J" e aplicar a função ''sample'' nele.
+<code rsplus>
+vetor <- rep(LETTERS[1:10])
 vetor
 sample(vetor)
-sample(vetor, replace=T)
+sample(vetor, replace = TRUE)
-sample(vetor,40,replace=T)
+sample(vetor, 40, replace = TRUE)
-sample(vetor,prob=c(0.1,0.2,0.05,0.05,0.2,0.1,0.05,0.05,0.1,0.1),replace=T)
+sample(vetor, prob = c(0.1,0.2,0.05,0.05,0.2,0.1,0.05,0.05,0.1,0.1),replace=TRUE)
-## O argumento ''prob'' é padronizado para somar um:
+sample(vetor, prob = c(1,2,0.5,0.5,2,1,0.5,0.5,1,1), replace = TRUE)
-sample(vetor,prob=c(1,2,0.5,0.5,2,1,0.5,0.5,1,1),replace=T) # o argumento prob pode somar mais que um
 </code>
-===== Dados de mandíbula de Chacal Dourado =====
-Vamos voltar aos dados de Chacal e à pergunta se há diferença no tamanho de mandíbulas entre machos e fêmeas
-<code>
+===== 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 tutorial [[02_tutoriais:tutorial6:start|]].
+<code rsplus>
 macho=c(120,107,110,116, 114, 111, 113,117,114,112)
 femea=c(110,111,107, 108,110,105,107,106,111,111)
@@ Linha 50: / 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.
-===== 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>
+<code rsplus>
 s1.mf=sample(mf)
 s1.mf
@@ Linha 67: / Linha 90: @@
 </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 criar um loop!!!!!
+<WRAP center round tip 60%>
+==== Ciclos de iteração ====
-<code>
+<code rsplus>
-result<-rep(NA,1000)
-result[1]<-diff(tapply(mf,sexo,mean))
+for(var in seq)
+{
+}
+</code>
+Onde:
+  * ''var'' é o nome sintético para uma variável
+  * ''seq'' vetor de valores que será assumido por ''var''
+  * ''{}'' expressões de procedimentos a serem repetidos
+</WRAP>
+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.
+<code rsplus>
+result <- rep(NA, 1000)
+result[1] <- diff(tapply(mf, sexo, mean))
 for(i in 2:1000)
 	{
-	dif.dados=diff(tapply(sample(mf),sexo,mean))
+	dif.dados <- diff(tapply(sample(mf), sexo, mean))
-	result[i]<-dif.dados
+	result[i] <- dif.dados
 	}
+</code>
+Um primeiro passo é fazer um gráfico com esse vetor de resultados:
+<code rsplus>
 hist(result)
-abline(v = result[1], col="red")
+abline(v = result[1], col = "red")
-abline(v = result[1]*-1, col="red")
+abline(v = result[1]*-1, col = "red")
 </code>
+==== Cálculo do p-valor ====
-===== Cálculo do P =====
+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.
-<code>
+Para o cálculo do p-valor da afirmação que há diferença temos:
-## Há diferença entre machos e fêmeas?
+<code rsplus>
 bicaudal=sum(result>=result[1]| result<=(result[1]*-1))
@@ Linha 97: / Linha 147: @@
 p.bi
-## Machos são maiores que as fêmeas?
+</code>
+Para se os machos são maiores que as fêmeas:
+<code rsplus>
 unicaudal=sum(result>=result[1])
 unicaudal
@@ Linha 104: / Linha 157: @@
 p.uni
 </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.
 ===== 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:
-<code>
+<code rsplus>
 macho
 macho.m
 </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:
-<code>
+Agora, fazemos um aleatorização deste vetor e calcular novamente a média:
+<code rsplus>
 mean(sample(macho))
 </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:
-<code>
+<code rsplus>
-smacho<-sample(macho, replace=TRUE)
+smacho <- sample(macho, replace = TRUE)
 mean(smacho)
-mean(sample(macho,replace=TRUE))
+mean(sample(macho, replace = TRUE))
-mean(sample(macho,replace=TRUE))
+mean(sample(macho, replace = TRUE))
 </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>
-nsim=100
+<code rsplus>
-resulta=rep(NA,nsim)
+nsim <- 100
+resulta <- rep(NA,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
 </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.
-<code>
+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>
 sort(resulta)
 sort(resulta)[6] ## o valore que deixa as 5 menores de fora
 sort(resulta)[95] ## o valore que deixa os 5 maiores de fora
 </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>
+<code rsplus>
 quantile(resulta, prob=c(0.05, 0.95))
 </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.
+/*
 ===== Função Vegas ======
@@ Linha 153: / Linha 235: @@
 Veja se você é capaz de entender o que a função faz a cada linha de comando e se estaria apto a explicá-la a outra pessoa.
-{{:bie5782:02_tutoriais:vegast.r|função vegas.t}}
+{{:02_tutoriais:vegast.r|função vegas.t}}
 Agora use a função para testar as hipóteses novamente!
+*/
 ===== 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 [[bie5782:03_apostila:08-simulacao#fnt__1|livro do Manly]] sobre permutação.
-A ideia era verificar se a ocorrência de taxa de tesourinhas (//Dermaptera//) estava mais correlacionada com a distribuição dos continentes atual ou antes da deriva continental.
+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 informação que partimos é do coeficiente de correlação da ocorrência de taxa 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.cef'':
+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.
-<code>
+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'':
+<code rsplus>
 data.coef<-matrix(c(NA, .30, .14, .23, .30, -0.04, 0.02, -0.09, NA, NA, .50,.50, .40, 0.04, 0.09, -0.06, NA, NA, NA, .54, .50, .11, .14, 0.05, rep(NA, 4), .61, .03,-.16, -.16, rep(NA, 5), .15, .11, .03, rep(NA, 6), .14, -.06, rep(NA, 7), 0.36, rep(NA, 8)), nrow=8, ncol=8)
 rownames(data.coef) <- c("Eur_Asia", "Africa", "Madag", "Orient", "Austr", "NewZea", "SoutAm", "NortAm")
@@ Linha 169: / Linha 253: @@
 </code>
-Foram usadas nesse estudo outras duas matrizes de distância, a primeira representando a distância atuais e a outra a distância geográfica antes da deriva continental das mesmas regiões biogeográficas.
+Foram usadas nesse estudo outras duas matrizes de distância, a primeira representando o número de eventos de dispersão de longas distâncias necessários para a conexão de populações na configuração atual dos continentes e a outra na configuração antes da deriva continental, entre as mesmas regiões biogeográficas.
-<code>
+<code rsplus>
 dist.atual<-matrix(c(NA,1,2,1,2,3,2,1, NA, NA, 1,2,3,4,3,2, NA, NA, NA,3,4,5,4,3, rep(NA, 4),1,2,3,2, rep(NA, 5), 1,4,3, rep(NA, 6), 5,4, rep(NA, 7), 1, rep(NA, 8)), nrow=8, ncol=8)
 dist.atual
@@ Linha 185: / Linha 269: @@
 </code>
-A primeira parte da análise dos dados e ver qual a correlação entre a matriz de correlação taxonômica e as distâncias geográficas (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 vamos calcular um coeficiente de correlação de //Pearson// entre as matrizes. Esse valor irá nos dizer se as duas matrizes estão correlacionadas, ou seja, os valores de uma variam na mesma direção da outra (+1), em direção contrária (-1) ou não são correlacionadas (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}}}$$
-<code>
+<code rsplus>
 cor12<-cor(as.vector(data.coef), as.vector(dist.atual), use="complete.obs")
 cor13<-cor(as.vector(data.coef), as.vector(dist.deriva), use="complete.obs")
@@ Linha 196: / Linha 280: @@
 </code>
-Ambos os valores de correlação estão nos dizendo que quanto maior a distância geográfica mais diferente é a composição de espécies de tesourinha. Além disso, que a correlação com as distâncias antes da deriva é mais forte. No caso, valores maiores em módulo já que a relação é de correlação negativa (aumento da distância diminui a similaridade florística).
+Ambos os valores de correlação estão nos dizendo que, quanto maior a distância geográfica mais diferente é a composição de espécies de tesourinha. Além disso, que a correlação com as distâncias antes da deriva é mais forte. No caso, valores maiores em módulo, já que a relação é de correlação negativa (aumento da distância diminui a similaridade florística).
-Agora precisamos calcular se esse valores de correlação poderiam ser atribuídos ao acaso. Para isso vamos fazer a permutação de uma das matrizes em e calcular o coeficientes de Pearson após essa permutação.
+Agora vamos calcular se esse valores de correlação poderiam ser atribuídos ao acaso. Para isso, vamos fazer a permutação de uma das matrizes e calcular o coeficientes de Pearson, após essa permutação.
-A permutação é simples, vamos mudar as colunas e linhas de lugares de maneira a aleatorizar os valores mas manter a estrutura subjacente ao dados. Uma maneira de fazer é:
+A permutação é simples, vamos mudar as colunas e linhas de lugares de maneira a aleatorizar os valores, mas manter a estrutura subjacente ao dados. Uma maneira de fazer é:
-<code>
+<code rsplus>
 data.sim<-data.coef # copia da matriz que será aleatorizada
 data.sim
@@ Linha 219: / Linha 303: @@
 cor12 ## correlação observada com a distancia atual
 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))
 str(res.cor)
@@ Linha 233: / Linha 320: @@
         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)
 par(mfrow=c(2,1))
@@ Linha 244: / Linha 338: @@
 p13=sum(res.cor[,2]<= res.cor[1,2])/(dim(res.cor)[1])
 p13
 </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 =====
+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:
+  * Manly B. F. J., 1997 Randomization, bootstrap and Monte Carlo methods in biology. 2nd Ed., Chapman and Hall, London
+  * GOTELLI, N. J. & G. R. GRAVES. 1996. Null models in ecology. Washington and London, Smithsonian Institution Press
+Caso tenham interesse pelo assunto sugiro iniciar por eles. O livro do Gotelli está esgotado, mas o autor disponibiliza o PDF em seu site.

ecoR

Ferramentas do usuário

Ferramentas do site

Diferenças

Ferramentas da página