Ferramentas do usuário
03_apostila:08-simulacao
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| Próxima revisão | Revisão anterior | ||
|
03_apostila:08-simulacao [2020/08/12 06:04] 127.0.0.1 edição externa |
03_apostila:08-simulacao [2020/10/04 20:49] adalardo |
||
|---|---|---|---|
| Linha 1: | Linha 1: | ||
| <WRAP tabs> | <WRAP tabs> | ||
| - | * [[bie5782:03_apostila:08-simulacao| Apostila]] | + | * [[03_apostila:08-simulacao| Apostila]] |
| - | * [[bie5782:02_tutoriais:tutorial9:start|Tutorial]] | + | * [[02_tutoriais:tutorial9:start|Tutorial]] |
| - | * [[bie5782:01_curso_atual:exerPermuta| Exercícios]] | + | * [[01_curso_atual:exerPermuta| Exercícios]] |
| </WRAP> | </WRAP> | ||
| ====== 8. Reamostragem e Simulação ====== | ====== 8. Reamostragem e Simulação ====== | ||
| - | O R é uma ferramenta poderosa para simular situações e compará-las com padrões obeservados. O assunto é extenso, e aqui apresentaremos apenas exemplos simples das modalidades mais usadas dem ecologia. Para quem quiser se se aprofundar, um bom começo é Manly (1997 ((Manly B. F. J., 1997 Randomization, bootstrap and Monte Carlo methods in biology. 2nd Ed., Chapman and Hall, London))). | + | O R é uma ferramenta poderosa para simular situações e compará-las com padrões observados. O assunto é extenso, e aqui apresentaremos apenas exemplos simples das modalidades mais usadas dem ecologia. Para quem quiser se se aprofundar, um bom começo é Manly (1997 ((Manly B. F. J., 1997 Randomization, bootstrap and Monte Carlo methods in biology. 2nd Ed., Chapman and Hall, London))). |
| ===== A Função ''sample'' ===== | ===== A Função ''sample'' ===== | ||
| Esta função reamostra os elementos de um vetor: | Esta função reamostra os elementos de um vetor: | ||
| - | <code> | + | <code rsplus> |
| > meu.vetor | > meu.vetor | ||
| [1] "A" "A" "B" "B" "C" "C" "D" "D" "E" "E" | [1] "A" "A" "B" "B" "C" "C" "D" "D" "E" "E" | ||
| Linha 35: | Linha 35: | ||
| Os dados que usaremos são numero de besouros no conteúdo estomacal de lagartos //Phrynosoma brevirostrae// machos e femeas (Manly,1997): | Os dados que usaremos são numero de besouros no conteúdo estomacal de lagartos //Phrynosoma brevirostrae// machos e femeas (Manly,1997): | ||
| - | <code> | ||
| - | <code> | + | <code rsplus> |
| > lagartos <- data.frame(sexo=factor(rep(c("m","f"),c(24,21))), | > lagartos <- data.frame(sexo=factor(rep(c("m","f"),c(24,21))), | ||
| + ncoleop=c(256,209,0,0,0,44,49,117,6,0,0,75,34,13,0,90,0,32,0,205,332,0,31,0, | + ncoleop=c(256,209,0,0,0,44,49,117,6,0,0,75,34,13,0,90,0,32,0,205,332,0,31,0, | ||
| Linha 45: | Linha 44: | ||
| Uma análise exploratória dos dados mostra que as médias amostrais são muito diferentes, mas as variâncias também, e a variação é muito grande, devido ao grande número de animais sem sinais de besouros no estômago: | Uma análise exploratória dos dados mostra que as médias amostrais são muito diferentes, mas as variâncias também, e a variação é muito grande, devido ao grande número de animais sem sinais de besouros no estômago: | ||
| - | <code> | + | <code rsplus> |
| > tapply(lagartos$ncoleop,lagartos$sexo,summary) | > tapply(lagartos$ncoleop,lagartos$sexo,summary) | ||
| $f | $f | ||
| Linha 73: | Linha 72: | ||
| Os gráficos de quantis téoricos indicam um forte desvio em relação à distribuição normal: | Os gráficos de quantis téoricos indicam um forte desvio em relação à distribuição normal: | ||
| - | <box center orange 70%>{{:bie5782:03_apostila:machos.jpg|}}</box> | + | <box center orange 70%>{{:03_apostila:machos.jpg|}}</box> |
| - | <box center orange 70%>{{:bie5782:03_apostila:femeas.jpg|}}</box> | + | <box center orange 70%>{{:03_apostila:femeas.jpg|}}</box> |
| Como os dados não atendem as premissas de um teste t, uma alternativa é substituí-lo por um teste de permutação. Como a hipótese nula é que as duas amostras vêm de populações com a mesma média, ela pode ser simulada permutando-se ao acaso os valores entre os sexos. | Como os dados não atendem as premissas de um teste t, uma alternativa é substituí-lo por um teste de permutação. Como a hipótese nula é que as duas amostras vêm de populações com a mesma média, ela pode ser simulada permutando-se ao acaso os valores entre os sexos. | ||
| Linha 82: | Linha 81: | ||
| A permutação é feita com a função ''sample'', que pode ser repetida com um //loop//, por meio da função ''for'': | A permutação é feita com a função ''sample'', que pode ser repetida com um //loop//, por meio da função ''for'': | ||
| - | <code> | + | <code rsplus> |
| > ##cria um vetor para armazenar os resultados | > ##cria um vetor para armazenar os resultados | ||
| > results <- c() | > results <- c() | ||
| Linha 93: | Linha 92: | ||
| Em apenas 15 das mil permutações a diferença absoluta das médias foi igual ou maior do que a observada: | Em apenas 15 das mil permutações a diferença absoluta das médias foi igual ou maior do que a observada: | ||
| - | <code> | + | <code rsplus> |
| > sum(results >= dif.obs) | > sum(results >= dif.obs) | ||
| [1] 15 | [1] 15 | ||
| Linha 99: | Linha 98: | ||
| Logo, a diferença observada é pouco provável sob a hipótese nula, o que nos permite rejeitá-la. Um gráfico de densidade probabilística dos valores das permutações ilustra isso: | Logo, a diferença observada é pouco provável sob a hipótese nula, o que nos permite rejeitá-la. Um gráfico de densidade probabilística dos valores das permutações ilustra isso: | ||
| - | <code> | + | <code rsplus> |
| > plot(density(results),xlab="Diferença Absoluta das Médias",ylab="Freq Relativa", main="") | > plot(density(results),xlab="Diferença Absoluta das Médias",ylab="Freq Relativa", main="") | ||
| > abline(v = dif.obs, col="red") | > abline(v = dif.obs, col="red") | ||
| Linha 105: | Linha 104: | ||
| </code> | </code> | ||
| - | <box orange center 100% >{{:bie5782:03_apostila:freq_relativa.jpg|}}</box> | + | <box orange center 100% >{{:03_apostila:freq_relativa.jpg|}}</box> |
| ==== Modelos Nulos em Ecologia de Comunidades ==== | ==== Modelos Nulos em Ecologia de Comunidades ==== | ||
| Linha 113: | Linha 112: | ||
| Vamos imaginar uma dessas matrizes, com seis espécies (linhas) e seis locais (colunas): | Vamos imaginar uma dessas matrizes, com seis espécies (linhas) e seis locais (colunas): | ||
| - | <code> | + | <code rsplus> |
| > cc2 | > cc2 | ||
| [,1] [,2] [,3] [,4] [,5] [,6] | [,1] [,2] [,3] [,4] [,5] [,6] | ||
| Linha 135: | Linha 134: | ||
| Para isso, basta permutar os valores dentro de cada linha: | Para isso, basta permutar os valores dentro de cada linha: | ||
| - | <code> | + | <code rsplus> |
| > apply(cc2,1,sample) | > apply(cc2,1,sample) | ||
| sp1 sp2 sp3 sp4 sp5 sp6 | sp1 sp2 sp3 sp4 sp5 sp6 | ||
| Linha 148: | Linha 147: | ||
| Ops! a função ''apply'' retorna sempre seus resultados em colunas, transpondo a matriz original. Para evitar isso, transpomos o resultado, voltando as espécies para as linhas: | Ops! a função ''apply'' retorna sempre seus resultados em colunas, transpondo a matriz original. Para evitar isso, transpomos o resultado, voltando as espécies para as linhas: | ||
| - | <code> | + | <code rsplus> |
| > t(apply(cc2,1,sample)) | > t(apply(cc2,1,sample)) | ||
| [,1] [,2] [,3] [,4] [,5] [,6] | [,1] [,2] [,3] [,4] [,5] [,6] | ||
| Linha 161: | Linha 160: | ||
| Agora basta repetir essa permutação muitas vezes e calcular o número máximo de espécies em co-ocorrência. Para esse cálculo, podemos criar uma função: | Agora basta repetir essa permutação muitas vezes e calcular o número máximo de espécies em co-ocorrência. Para esse cálculo, podemos criar uma função: | ||
| - | <code> | + | <code rsplus> |
| > ##Uma funcao para calcular o maximo de especies em co-ocorrencia | > ##Uma funcao para calcular o maximo de especies em co-ocorrencia | ||
| > max.c <- function(x){ max(apply(x,2,sum)) } | > max.c <- function(x){ max(apply(x,2,sum)) } | ||
| Linha 171: | Linha 170: | ||
| E então repetimos as permutações 1000 vezes, inserindo-as em um //loop//: | E então repetimos as permutações 1000 vezes, inserindo-as em um //loop//: | ||
| - | <code> | + | <code rsplus> |
| ##Um vetor para guardar os resultados | ##Um vetor para guardar os resultados | ||
| > results <- c() | > results <- c() | ||
| Linha 182: | Linha 181: | ||
| O número de randomizações que teve um máximo de espécies por local menor ou igual ao observado foi alto (21,6%), portanto a observação de que no máximo três espécies ocorrem no mesmo local não é evidência de que haja um limite a esse número: | O número de randomizações que teve um máximo de espécies por local menor ou igual ao observado foi alto (21,6%), portanto a observação de que no máximo três espécies ocorrem no mesmo local não é evidência de que haja um limite a esse número: | ||
| - | <code> | + | <code rsplus> |
| > sum(result.c<=max.coc) | > sum(result.c<=max.coc) | ||
| [1] 216 | [1] 216 | ||
| Linha 198: | Linha 197: | ||
| A aplicação mais simples deste princípio é estimar intervalos de confiança por reamostragem. Para exemplificar, usaremos o conjunto de dados BCI, um dataframe de abundâncias de espécies de árvores na parcela permanente de Barro Colorado nas (linhas), por subparcelas de 1 ha: | A aplicação mais simples deste princípio é estimar intervalos de confiança por reamostragem. Para exemplificar, usaremos o conjunto de dados BCI, um dataframe de abundâncias de espécies de árvores na parcela permanente de Barro Colorado nas (linhas), por subparcelas de 1 ha: | ||
| - | <code> | + | <code rsplus> |
| ##Esse objeto de dados está no pacote "vegan" | ##Esse objeto de dados está no pacote "vegan" | ||
| > data(BCI, package="vegan") | > data(BCI, package="vegan") | ||
| Linha 210: | Linha 209: | ||
| Vamos criar uma função para calcular o índice de Shannon a partir de um vetor de abundâncias de espécies, e aplicar a cada parcela: | Vamos criar uma função para calcular o índice de Shannon a partir de um vetor de abundâncias de espécies, e aplicar a cada parcela: | ||
| - | <code> | + | <code rsplus> |
| > H <- function(x){ | > H <- function(x){ | ||
| + y <- x[x>0] | + y <- x[x>0] | ||
| Linha 220: | Linha 219: | ||
| O vetor H.BCI contém os índices de Shannon de cada parcela de 1 ha. Sua distribuição não parece ser normal: | O vetor H.BCI contém os índices de Shannon de cada parcela de 1 ha. Sua distribuição não parece ser normal: | ||
| - | <code> | + | <code rsplus> |
| > hist(H.BCI,xlab="Índice de Shannon") | > hist(H.BCI,xlab="Índice de Shannon") | ||
| > qqnorm(H.BCI);qqline(H.BCI) | > qqnorm(H.BCI);qqline(H.BCI) | ||
| </code> | </code> | ||
| - | <box center orange 100%>{{:bie5782:03_apostila:histogram_bci.jpg|}}{{:bie5782:03_apostila:normal_qq_plot.jpg|}}</box> | + | <box center orange 100%>{{:03_apostila:histogram_bci.jpg|}}{{:03_apostila:normal_qq_plot.jpg|}}</box> |
| Vamos simular uma amostra deste conjuntos de parcelas, e estimar o intervalo de confiança da diversidade média por parcela: | Vamos simular uma amostra deste conjuntos de parcelas, e estimar o intervalo de confiança da diversidade média por parcela: | ||
| - | <code> | + | <code rsplus> |
| > ##Uma amostra ao acaso de 20 parcelas | > ##Uma amostra ao acaso de 20 parcelas | ||
| > bci.sample <- BCI[sample(1:nrow(BCI),20),] | > bci.sample <- BCI[sample(1:nrow(BCI),20),] | ||
| Linha 252: | Linha 251: | ||
| No caso, esse intervalo não inclui o zero, o que rejeita a hipótese nula. Mas o que interessa aqui é que o intervalo inclui a média da população de parcelas de onde veio a amostra, que é : | No caso, esse intervalo não inclui o zero, o que rejeita a hipótese nula. Mas o que interessa aqui é que o intervalo inclui a média da população de parcelas de onde veio a amostra, que é : | ||
| - | <code> | + | <code rsplus> |
| > mean(H.BCI) | > mean(H.BCI) | ||
| [1] 3.82084 | [1] 3.82084 | ||
| Linha 261: | Linha 260: | ||
| Como a premissa de normalidade parece não ser verdadeira, uma solução é calcular um intervalo por //bootstrap//. Para isso, basta re-amostrar com reposição a amostra, e calcular a estatística de interesse, no caso a média dos índices de Shannon de cada parcela: | Como a premissa de normalidade parece não ser verdadeira, uma solução é calcular um intervalo por //bootstrap//. Para isso, basta re-amostrar com reposição a amostra, e calcular a estatística de interesse, no caso a média dos índices de Shannon de cada parcela: | ||
| - | <code> | + | <code rsplus> |
| ##Um vetor para armazenar os resultados | ##Um vetor para armazenar os resultados | ||
| > results <- c() | > results <- c() | ||
| Linha 269: | Linha 268: | ||
| A diferença entre a média da amostra original e a média das aleatorizações estima o viés, no caso muito pequeno: | A diferença entre a média da amostra original e a média das aleatorizações estima o viés, no caso muito pequeno: | ||
| - | <code> | + | <code rsplus> |
| > mean(results) | > mean(results) | ||
| [1] 3.853739 | [1] 3.853739 | ||
| Linha 278: | Linha 277: | ||
| O intervalo de confiança por quantis a 95% está compreendida entre o valor abaixo do qual há apenas 2,5% das randomizaçôes e o valor acima do qual restam apenas 2,5% das randomizações. Portanto, são os quantis a 2,5% e 97,5% dos mil valores simulados. Usamos a função ''quantile'' para obtê-los: | O intervalo de confiança por quantis a 95% está compreendida entre o valor abaixo do qual há apenas 2,5% das randomizaçôes e o valor acima do qual restam apenas 2,5% das randomizações. Portanto, são os quantis a 2,5% e 97,5% dos mil valores simulados. Usamos a função ''quantile'' para obtê-los: | ||
| - | <code> | + | <code rsplus> |
| > quantile(results, probs=c(0.025,0.975)) | > quantile(results, probs=c(0.025,0.975)) | ||
| 2.5% 97.5% | 2.5% 97.5% | ||
| Linha 293: | Linha 292: | ||
| Para ilustrar a forma dessas distribuições, podemos fazer o gráfico dos valores teóricos de probabilidade, no intervalo de zero a 25: | Para ilustrar a forma dessas distribuições, podemos fazer o gráfico dos valores teóricos de probabilidade, no intervalo de zero a 25: | ||
| - | <code> | + | <code rsplus> |
| > x <- 0:25 | > x <- 0:25 | ||
| > y1 <- dnbinom(x,mu=2, size=2) | > y1 <- dnbinom(x,mu=2, size=2) | ||
| Linha 304: | Linha 303: | ||
| </code> | </code> | ||
| - | <box orange center 100%>{{:bie5782:03_apostila:probabilidadexvalores.png|}}</box> | + | <box orange center 100%>{{:03_apostila:probabilidadexvalores.png|}}</box> |
| As premissas do teste t são que as amostras vêm de populações com distribuições normais que têm variâncias iguais (pelo menos aproximadamente). A distribuição binomial negativa tem uma forma muito diferente da normal, e nela a variância é uma função da média. Assim, duas distribuições binomiais negativas com médias diferentes, como as ilustradas acima, terão variâncias diferentes. Além disso, a distribuição binomial negativa é **discreta**, enquanto a normal é **contínua**. | As premissas do teste t são que as amostras vêm de populações com distribuições normais que têm variâncias iguais (pelo menos aproximadamente). A distribuição binomial negativa tem uma forma muito diferente da normal, e nela a variância é uma função da média. Assim, duas distribuições binomiais negativas com médias diferentes, como as ilustradas acima, terão variâncias diferentes. Além disso, a distribuição binomial negativa é **discreta**, enquanto a normal é **contínua**. | ||
| Linha 312: | Linha 311: | ||
| Para simular a comparação de duas amostras de tamanho dez, usando um //loop// temos: | Para simular a comparação de duas amostras de tamanho dez, usando um //loop// temos: | ||
| - | <code> | + | <code rsplus> |
| > ##Um vetor fator com 20 elementos, para representar cada elemento das duas amostras | > ##Um vetor fator com 20 elementos, para representar cada elemento das duas amostras | ||
| > tratamento <- factor(x=rep(c("a","b"),each=10)) | > tratamento <- factor(x=rep(c("a","b"),each=10)) | ||
| Linha 330: | Linha 329: | ||
| Vamos agora avaliar o que acontece se tomamos amostras de populações com distribuição normal, mantendo as mesmas variâncias: | Vamos agora avaliar o que acontece se tomamos amostras de populações com distribuição normal, mantendo as mesmas variâncias: | ||
| - | <code> | + | <code rsplus> |
| > for(i in 1:10000){ | > for(i in 1:10000){ | ||
| + variaveis <- c(rnorm(n=10,mean=2,sd=2), | + variaveis <- c(rnorm(n=10,mean=2,sd=2), | ||
| Linha 342: | Linha 341: | ||
| Como poderíamos esperar, a força estimada foi maior, pois uma das premissas do teste agora é satisfeita. No entanto, a premissa de igualdade de variâncias ainda não é cumprida. Se simulamos amostras com variâncias iguais à média das duas variâncias usadas nas simulações anteriores, temos um ganho na força do teste, embora pequeno: | Como poderíamos esperar, a força estimada foi maior, pois uma das premissas do teste agora é satisfeita. No entanto, a premissa de igualdade de variâncias ainda não é cumprida. Se simulamos amostras com variâncias iguais à média das duas variâncias usadas nas simulações anteriores, temos um ganho na força do teste, embora pequeno: | ||
| - | <code> | + | <code rsplus> |
| > ## Desvio-padrão médio correspondente às variâncias de 7,5 e 4 | > ## Desvio-padrão médio correspondente às variâncias de 7,5 e 4 | ||
| > sd1 <- mean(sqrt(c(7.5,4))) | > sd1 <- mean(sqrt(c(7.5,4))) | ||
| Linha 364: | Linha 363: | ||
| <box left 100% | //**Exercício:** Que simulação é esta?//> | <box left 100% | //**Exercício:** Que simulação é esta?//> | ||
| Descubra que simulação é feita com código abaixo: | Descubra que simulação é feita com código abaixo: | ||
| - | <code> | + | <code rsplus> |
| > m1 <- matrix(data=rnbinom(n=10000,size=2,mu=2),nrow=10,ncol=1000) | > m1 <- matrix(data=rnbinom(n=10000,size=2,mu=2),nrow=10,ncol=1000) | ||
03_apostila/08-simulacao.txt · Última modificação: 2020/10/04 20:49 por adalardo
Exceto onde for informado ao contrário, o conteúdo neste wiki está sob a seguinte licença: CC Attribution-Noncommercial-Share Alike 4.0 International
