Ferramentas do usuário

Ferramentas do site


02_tutoriais:tutorial9:start

Diferenças

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

Link para esta página de comparações

Ambos lados da revisão anterior Revisão anterior
Próxima revisão
Revisão anterior
02_tutoriais:tutorial9:start [2022/06/21 17:21]
adalardo [Criando ciclos de eventos]
02_tutoriais:tutorial9:start [2024/09/11 12:48] (atual)
Linha 14: Linha 14:
 Especificamente para o testes de hipótese por aleatorização ou randomização dos dados, temos ao menos as seguintes etapas:  ​ 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%> <WRAP center round box 60%>
-  - Definir a estatística de interesse (EI)+  - Definir a estatística de interesse (**EI**)
   - Calcular a estatística de interesse a partir dos dados   - Calcular a estatística de interesse a partir dos dados
   - Estabelecer o cenário nulo   - Estabelecer o cenário nulo
   - Simular cenário nulo   - Simular cenário nulo
-  - Calcular a EI no cenário nulo (pseudovalor)+  - Calcular a **EI** no cenário nulo (pseudovalor)
   - Produzir a distribuição dos pseudovalores   - Produzir a distribuição dos pseudovalores
-  - Posicionar a EI observada na distribuição dos pseudovalores+  - Posicionar a **EI** observada na distribuição dos pseudovalores
   - Calcular o p-valor   - Calcular o p-valor
 </​WRAP>​ </​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()''​. ​ 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>
-macho=c(120,​107,​110,​116,​ 114, 111, 113,​117,​114,​112) +macho <- c(120,​107,​110,​116,​ 114, 111, 113,​117,​114,​112) 
-femea=c(110,​111,​107,​ 108,​110,​105,​107,​106,​111,​111)+femea <- c(110,​111,​107,​ 108,​110,​105,​107,​106,​111,​111)
 macho macho
 femea femea
-sexo=rep(c("​macho",​ "​femea"​),​ each=10)+sexo <- rep(c("​macho",​ "​femea"​),​ each=10)
 sexo  sexo 
-mf=c(macho,​femea)+mf <- c(macho,​femea)
 mf mf
-macho.m=mean(macho)+macho.m ​<- mean(macho)
 macho.m macho.m
-femea.m=mean(femea)+femea.m ​<- mean(femea)
 femea.m femea.m
-macho.m-femea.m +macho.m - femea.m 
-dif.mf=diff(tapply(mf,​sexo,​mean))+dif.mf ​<- diff(tapply(mf,​sexo,​mean))
 dif.mf dif.mf
 </​code>​ </​code>​
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)
 s1.mf s1.mf
 diff(tapply(s1.mf,​sexo,​mean)) diff(tapply(s1.mf,​sexo,​mean))
 ##+1 ##+1
-s2.mf=sample(mf)+s2.mf <- sample(mf)
 s2.mf s2.mf
 diff(tapply(s2.mf,​sexo,​mean)) diff(tapply(s2.mf,​sexo,​mean))
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.
  
  
 <code rsplus> <code rsplus>
-result<​-rep(NA,​1000) +result <- rep(NA, 1000) 
-result[1]<​-diff(tapply(mf,​sexo,​mean))+result[1] <- diff(tapply(mf,​ sexo, mean))
 for(i in 2:1000) 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) 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>​ </​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. ​  
 + 
 +Para o cálculo ​do p-valor da afirmação que há diferença temos:
  
 <code rsplus> <code rsplus>
-## Há diferença entre machos e fêmeas? 
  
-bicaudal=sum(result>​=result[1]| result<​=(result[1]*-1))+bicaudal ​<- sum(result >= result[1]| result <= (result[1]*-1))
 bicaudal bicaudal
 length(result) length(result)
-p.bi=bicaudal/​length(result)+p.bi <- bicaudal/​length(result)
 p.bi p.bi
  
-## Machos são maiores que as fêmeas?+</​code>​
  
-unicaudal=sum(result>​=result[1])+Para a hipótese direcionada se os machos tem mandíbulas maiores que as fêmeas: 
 + 
 +<code rsplus>​ 
 +unicaudal ​<- sum(result >= result[1])
 unicaudal unicaudal
-p.uni=unicaudal/​length(result)+p.uni <- unicaudal/​length(result)
 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. ​
 +
 +
  
 ===== Bootstrap ===== ===== Bootstrap =====
  
-Vamos agora pegar o mesmo exemplo anterior ​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 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 155: 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 duas últimas linhas de comando produzem valores diferentesapesar de serem identicas. 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.as mesmas 
-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 portantopodemos 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 ​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%...). ​Para calcular um intervalo de 90% podemos ordenar o vetor verificar ​quais valores estão nos limite das posições ​com 5% de cada lado. Como o vetor tem cem posições só precisamos olhar a posição ''​6''​ e a posição ''​95''​ 
 +  
 <code rsplus> <code rsplus>
 sort(resulta) sort(resulta)
 sort(resulta)[6] ## o valore que deixa as 5 menores de fora sort(resulta)[6] ## o valore que deixa as 5 menores de fora
 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>​ 
-Podemos também usar a função //​quantile()//​ definindo os quantis de interesse: 
-<code rsplus> 
-quantile(resulta,​ prob=c(0.05,​ 0.95)) 
 </​code>​ </​code>​
  
  
-===== Simula T ===== +Podemos ​também ​usar a função ​''​quantile()''​ definindo os quantis ​de interesse:
- +
-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> <code rsplus>
-x11(width = 10height ​10) +quantile(resultaprob c(0.05, 0.95))
-source("​simulaT.r") +
-simulaT(macho,​ femea, teste = "​maior",​ anima = TRUE)+
 </​code>​ </​code>​
  
-As funções ​não precisam ser consideradas procedimento abstrato no qual não temos acessoUma função similar ​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çãoNa próxima aula iremos entender como incorporar um procedimento em uma função. O primeiro passo é saber executar ​procedimento em linhas de código, como fizemos no início do tutorial+Definitivamente,​ fazer só 100 simulações, ​não parece adequadoExistem muitos arranjos possíveis de 10 elementos reamostrados com reposição((combinatória com reposição ​de 10 elementos ​em conjuntos ​de 10 é da ordem de 92378))Refaça ​o código ​com 1000 (mil) iterações e recalcule o intervalo.
  
  
Linha 218: Linha 247:
  
 ===== Tesourinha e a deriva continental ===== ===== Tesourinha e a deriva continental =====
-Vamos agora reproduzir ​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'':​
 <code rsplus> <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)+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"​) rownames(data.coef) <- c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​)
 colnames(data.coef) <- c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​) colnames(data.coef) <- c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​)
Linha 231: Linha 261:
  
 <code rsplus> <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 <- 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 dist.atual
-dist.deriva<​- matrix(c(NA,​1,​2,​1,​2,​3,​2,​1,​ NA, NA, 1,​1,​1,​2,​1,​2,​ NA, NA, NA,​1,​1,​2,​2,​3,​ rep(NA, 4),1,2,2,2, rep(NA, 5), 1,2,3, rep(NA, 6), 3,4, rep(NA, 7), 1,  rep(NA, 8)), nrow=8, ncol=8)+dist.deriva <- matrix(c(NA,​1,​2,​1,​2,​3,​2,​1,​ NA, NA, 1,​1,​1,​2,​1,​2,​ NA, NA, NA,​1,​1,​2,​2,​3,​ rep(NA, 4),1,2,2,2, rep(NA, 5), 1,2,3, rep(NA, 6), 3,4, rep(NA, 7), 1,  rep(NA, 8)), nrow=8, ncol=8)
 # colocando nomes nas matrizes # colocando nomes nas matrizes
 rownames(dist.atual) <- colnames(dist.atual)<​- c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​) rownames(dist.atual) <- colnames(dist.atual)<​- c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​)
  
-colnames(dist.deriva)<​- rownames(dist.deriva)<​- ​ c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​)+colnames(dist.deriva) <- rownames(dist.deriva)<​- ​ c("​Eur_Asia",​ "​Africa",​ "​Madag",​ "​Orient",​ "​Austr",​ "​NewZea",​ "​SoutAm",​ "​NortAm"​)
 # olhando as matrizes # olhando as matrizes
 dist.atual dist.atual
Linha 245: Linha 275:
  
 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}}}$$
 +
 <code rsplus> <code rsplus>
-cor12<​-cor(as.vector(data.coef),​ as.vector(dist.atual),​ use="​complete.obs"​) +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"​)+cor13 <- cor(as.vector(data.coef),​ as.vector(dist.deriva),​ use="​complete.obs"​)
 cor12 ## correlação com a distancia atual cor12 ## correlação com a distancia atual
 cor13 ## correlação com a distancia antes da deriva cor13 ## correlação com a distancia antes da deriva
Linha 261: Linha 292:
  
 <code rsplus> <code rsplus>
-data.sim<​-data.coef # copia da matriz que será aleatorizada+data.sim <- data.coef # copia da matriz que será aleatorizada
 data.sim ​ data.sim ​
  
Linha 269: Linha 300:
 data.sim # olhando a matriz data.sim # olhando a matriz
 data.sim[8:​1,​ 8:1] # uma matriz baguncada mas que mantem certa estrutura data.sim[8:​1,​ 8:1] # uma matriz baguncada mas que mantem certa estrutura
-sim.pos<​-sample(1:​8) # posicoes permutadas ​+sim.pos <- sample(1:8) # posicoes permutadas ​
 sim.pos sim.pos
-data.sim<​-data.sim[sim.pos,​ sim.pos] # aqui uma matriz verdadeiramente permutada +data.sim <- data.sim[sim.pos,​ sim.pos] # aqui uma matriz verdadeiramente permutada 
-cor12.sim<​-cor(as.vector(data.sim),​ as.vector(dist.atual),​ use="​pairwise.complete.obs"​) +cor12.sim <- cor(as.vector(data.sim),​ as.vector(dist.atual),​ use="​pairwise.complete.obs"​) 
-cor13.sim<​-cor(as.vector(data.sim),​ as.vector(dist.deriva),​ use="​pairwise.complete.obs"​)+cor13.sim <- cor(as.vector(data.sim),​ as.vector(dist.deriva),​ use="​pairwise.complete.obs"​)
 cor12.sim ​ cor12.sim ​
 cor13.sim cor13.sim
 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.  
-res.cor=data.frame(sim12=rep(NA,​ 5000), sim13=rep(NA,​5000))+ 
 +<code rsplus> 
 + 
 +res.cor ​<- data.frame(sim12=rep(NA,​ 5000), sim13=rep(NA,​5000))
 str(res.cor) str(res.cor)
-res.cor[1,​]<​-c(cor12,​ cor13)+res.cor[1,] <- c(cor12, cor13)
 str(res.cor) str(res.cor)
 for(s in 2:5000) for(s in 2:5000)
     {     {
-        sim.pos<​-sample(1:​8) +        sim.pos <- sample(1:​8) 
-        data.sim<​-data.sim[sim.pos,​ sim.pos] +        data.sim <- data.sim[sim.pos,​ sim.pos] 
-        res.cor[s,​1]<​-cor(as.vector(data.sim),​ as.vector(dist.atual),​ use="​pairwise.complete.obs"​) +        res.cor[s,​1] <- cor(as.vector(data.sim),​ as.vector(dist.atual),​ use="​pairwise.complete.obs"​) 
-        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))
 hist(res.cor[,​1]) hist(res.cor[,​1])
-abline(v=res.cor[1,​1],​ col="​red"​)+abline(v = res.cor[1,​1],​ col = "​red"​)
 hist(res.cor[,​2]) hist(res.cor[,​2])
-abline(v=res.cor[1,​2],​ col="​red"​)+abline(v = res.cor[1,​2],​ col = "​red"​)
 #### calculando o P ########### #### calculando o P ###########
-p12=sum(res.cor[,​1]<​= res.cor[1,​1])/​(dim(res.cor)[1])+p12 <- sum(res.cor[,​1] <= res.cor[1,​1])/​(dim(res.cor)[1])
 p12 p12
-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.1655842890.txt.gz · Última modificação: 2022/06/21 17:21 por adalardo