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: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 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 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 diferentes, apesar 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 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%...). Para calcular um intervalo de 90% podemos ordenar o vetor e 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 = 10, height = 10) | + | quantile(resulta, prob = 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 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. | + | Definitivamente, fazer só 100 simulações, não parece adequado. Existem muitos arranjos possíveis de 10 elementos reamostrados com reposição((a 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 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'': | ||
<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: |