Diferenças

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

--- 02_tutoriais:tutorial9:start [2020/09/23 17:06]
adalardo code
+++ 02_tutoriais:tutorial9:start [2023/09/12 10:46]
@@ Linha 1: / Linha 1: @@
-<WRAP tabs>
-  * [[bie5782:02_tutoriais:tutorial9:start|Tutorial]]
-  * [[bie5782:01_curso_atual:exerPermuta| Exercícios]]
-  * [[bie5782: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 rsplus>
-vetor=rep(LETTERS[1:10])
-vetor
-sample(vetor)
-sample(vetor, replace=T)
-sample(vetor,40,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=T)
-## 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=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 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)
-macho
-femea
-sexo=rep(c("macho", "femea"), each=10)
-sexo
-mf=c(macho,femea)
-mf
-macho.m=mean(macho)
-macho.m
-femea.m=mean(femea)
-femea.m
-macho.m-femea.m
-dif.mf=diff(tapply(mf,sexo,mean))
-dif.mf
-</code>
-**PERGUNTAS:**
-  * Essa diferença entre as médias é significativa?
-  * Qual minha incerteza ao afirmar que essas médias são diferentes?
- 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 =====
-<code rsplus>
-s1.mf=sample(mf)
-s1.mf
-diff(tapply(s1.mf,sexo,mean))
-##+1
-s2.mf=sample(mf)
-s2.mf
-diff(tapply(s2.mf,sexo,mean))
-##+2
-diff(tapply(sample(mf),sexo,mean))
-##+3
-diff(tapply(sample(mf),sexo,mean))
-##+1000
-### e agora? fazer na mão 1000 vezes? ###
-</code>
-===== Criando ciclos de eventos =====
-Vamos criar um loop!!!!!
-<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))
-	result[i]<-dif.dados
-	}
-hist(result)
-abline(v = result[1], col="red")
-abline(v = result[1]*-1, col="red")
-</code>
-===== Cálculo do P =====
-<code rsplus>
-## Há diferença entre machos e fêmeas?
-bicaudal=sum(result>=result[1]| result<=(result[1]*-1))
-bicaudal
-length(result)
-p.bi=bicaudal/length(result)
-p.bi
-## Machos são maiores que as fêmeas?
-unicaudal=sum(result>=result[1])
-unicaudal
-p.uni=unicaudal/length(result)
-p.uni
-</code>
-===== Bootstrap =====
-Vamos agora pegar o mesmo exemplo anterior e estimar o intervalo de confiança da média dos machos do chacal dourado.
-Primeiro vamos ver novamente esses dados e  sua média:
-<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 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 rsplus>
-smacho<-sample(macho, replace=TRUE)
-mean(smacho)
-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.
-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.
-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>
-nsim=100
-resulta=rep(NA,nsim)
-for(i in 1:nsim)
- {
-  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 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:
-<code rsplus>
-quantile(resulta, prob=c(0.05, 0.95))
-</code>
-===== Função Vegas ======
-Em aula, se houve tempo, construímos uma função que  automatiza a sequência de comandos da primeira parte desse tutorial onde testamos as hipóteses: (1) da mandíbulas de chacais machos e fêmeas serem diferentes e (2) a mandíbula de machos serem maiores que a das fêmeas, em média.
-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}}
-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.
-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'':
-<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")
-colnames(data.coef) <- c("Eur_Asia", "Africa", "Madag", "Orient", "Austr", "NewZea", "SoutAm", "NortAm")
-data.coef
-</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.
-<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
-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
-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")
-# olhando as matrizes
-dist.atual
-dist.deriva
-</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).
-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).
- $$ 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>
-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")
-cor12 ## correlação com a distancia atual
-cor13 ## correlação com a distancia antes da deriva
-</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).
-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.
-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 rsplus>
-data.sim<-data.coef # copia da matriz que será aleatorizada
-data.sim
-# preenchendo o triangulo superior da matriz com os dados correspondentes do triangulo inferior
-data.sim[upper.tri(data.sim)] <- t(data.coef)[(upper.tri(data.coef))]
-data.sim # olhando a matriz
-data.sim[8:1, 8:1] # uma matriz baguncada mas que mantem certa estrutura
-sim.pos<-sample(1:8) # posicoes permutadas
-sim.pos
-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")
-cor13.sim<-cor(as.vector(data.sim), as.vector(dist.deriva), use="pairwise.complete.obs")
-cor12.sim
-cor13.sim
-cor12 ## correlação observada com a distancia atual
-cor13 ## correlação observada com a distancia antes da deriva
-########################################################
-### Repetir a simulação muitas vezes ###################
-#######################################################
-res.cor=data.frame(sim12=rep(NA, 5000), sim13=rep(NA,5000))
-str(res.cor)
-res.cor[1,]<-c(cor12, cor13)
-str(res.cor)
-for(s in 2:5000)
-    {
-        sim.pos<-sample(1:8)
-        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,2]<-cor(as.vector(data.sim), as.vector(dist.deriva), use="pairwise.complete.obs")
-    }
-str(res.cor)
-par(mfrow=c(2,1))
-hist(res.cor[,1])
-abline(v=res.cor[1,1], col="red")
-hist(res.cor[,2])
-abline(v=res.cor[1,2], col="red")
-#### calculando o P ###########
-p12=sum(res.cor[,1]<= res.cor[1,1])/(dim(res.cor)[1])
-p12
-p13=sum(res.cor[,2]<= res.cor[1,2])/(dim(res.cor)[1])
-p13
-</code>

ecoR

Ferramentas do usuário

Ferramentas do site

Diferenças

Ferramentas da página