* [[bie5782:02_tutoriais:tutorial9:start|Tutorial]] * [[bie5782:01_curso_atual:exerPermuta| Exercícios]] * [[bie5782:03_apostila:08-simulacao| Apostila]] ====== 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. 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 ===== 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 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 **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 ===== 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? ### ===== Criando ciclos de eventos ===== Vamos criar um loop!!!!! 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") ===== Cálculo do P ===== ## 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 ===== 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: macho macho.m 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: mean(sample(macho)) 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: smacho<-sample(macho, replace=TRUE) mean(smacho) mean(sample(macho,replace=TRUE)) mean(sample(macho,replace=TRUE)) 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: nsim=100 resulta=rep(NA,nsim) for(i in 1:nsim) { resulta[i]<-mean(sample(macho, replace=TRUE)) } ## veja os valores calculados resulta 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. 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 Podemos também usar a função //quantile()// definindo os quantis de interesse: quantile(resulta, prob=c(0.05, 0.95)) ===== 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'': 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 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. 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 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}}}$$ 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 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 é: 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