8. Reamostragem e Simulação [ecoR]

Essa página está em modo somente de leitura. Você pode visualizar a fonte, mas não alterá-la. Informe-se com o administrador do Wiki, caso você ache que isso está incorreto.
<WRAP tabs>
  * [[03_apostila:08-simulacao| Apostila]]
  * [[02_tutoriais:tutorial9:start|Tutorial]]
  * [[01_curso_atual:exerPermuta| Exercícios]] 

</WRAP>
====== 8. Reamostragem e Simulação ======

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'' =====
Esta função reamostra os elementos de um vetor:
<code rsplus>
> meu.vetor
 [1] "A" "A" "B" "B" "C" "C" "D" "D" "E" "E"
> meu.vetor.bagunçado <- sample(meu.vetor)
> meu.vetor.bagunçado
 [1] "C" "A" "A" "D" "B" "E" "E" "C" "D" "B"
</code>

Por //default// a função retorna um novo vetor de mesmo comprimento, com os elementos permutados ao acaso. Portanto, é uma amostragem **sem reposição**. A função tem argumentos que definem se a amostragem é ou não com reposição, o tamanho do vetor resultante, entre outros. Consulte a página de ajuda.

==== Exercícios ====

<box left | //**Exercício:** Quem Precisa da Caixa Econômica?//>
Você pode simular sorteios das loterias com a função sample. Consulte a ajuda da função para descobrir como simular um sorteio como o feito para o jogo da sena.
</box>

===== Reamostragens sem Reposição: Permutações =====
Nesse procedimentos apenas "embaralhamos" os elementos de um ou mais vetores, o que podemos fazer gerando um vetor de mesmo tamanho, cujos elementos são amostrados ao acaso **sem reposição**. Esse procedimento é chamado de permutação ao acaso.


==== Testes de Permutação ====
Uma das aplicações simples das permutações ao acaso são testes de comparação de médias entre grupos, como o teste t, como mostramos abaixo.

Os dados que usaremos são numero de besouros no conteúdo estomacal de lagartos //Phrynosoma brevirostrae// machos e femeas (Manly,1997):

<code rsplus>
> 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,
+                           0,89,0,0,0,163,289,3,843,0,158,443,311,232,179,179,19,142,100,0,432))
</code>

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 rsplus>
> tapply(lagartos$ncoleop,lagartos$sexo,summary)
$f
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    0.0     0.0   142.0   170.6   232.0   843.0

$m
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   0.00    0.00   22.00   62.21   78.75  332.00

##Médias
> medias <- tapply(lagartos$ncoleop,lagartos$sexo,mean)
> medias
        f         m
170.57143  62.20833
##Diferença absoluta entre as médias
> dif.obs <- abs( diff( as.vector(medias) ) )
> dif.obs
[1] 108.3631

##Variâncias
> tapply(lagartos$ncoleop,lagartos$sexo,var)
        f         m
43533.557  8855.911
</code>

Os gráficos de quantis téoricos indicam um forte desvio em relação à distribuição normal:

<box center orange 70%>{{:03_apostila:machos.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. 

Calculamos então um índice de diferenças entre as médias das amostras que deve ser similar ao obtido, caso a hipótese nula esteja correta. Repetindo esse procedimento milhares de vezes, podemos estimar a chance de um valor igual ou maior que o observado ter ocorrido mesmo se as duas amostras vêm de populações com médias diferentes.

A permutação é feita com a função ''sample'', que pode ser repetida com um //loop//, por meio da função ''for'':

<code rsplus>
> ##cria um vetor para armazenar os resultados
> results <- c()
> ##Permuta os valores das medidas, calcula a diferença absoluta entre as médias e 
> ##armazena no vetor "results". Repete a operação mil vezes
> for(i in 1:1000){
+ results[i] <- abs( diff( tapply(sample(lagartos$ncoleop),lagartos$sexo,mean) ) )
+ }
</code>

Em apenas 15 das mil permutações a diferença absoluta das médias foi igual ou maior do que a observada:
<code rsplus>
> sum(results >= dif.obs)
[1] 15
</code>

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 rsplus>
> plot(density(results),xlab="Diferença Absoluta das Médias",ylab="Freq Relativa", main="")
> abline(v = dif.obs, col="red")

</code>

<box orange center 100% >{{:03_apostila:freq_relativa.jpg|}}</box>

==== Modelos Nulos em Ecologia de Comunidades ====

Modelos nulos em ecologia de comunidades ((ver GOTELLI, N. J. & G. R. GRAVES. 1996. Null models in ecology. Washington and London, Smithsonian Institution Press)) muitas vezes envolvem permutação de elementos de uma matriz, em geral de ocorrência ou abundância de espécies por local. 

Vamos imaginar uma dessas matrizes, com seis espécies (linhas) e seis locais (colunas):

<code rsplus>
> cc2
    [,1] [,2] [,3] [,4] [,5] [,6]
sp1    1    1    0    0    0    0
sp2    1    1    1    0    0    0
sp3    0    0    1    1    0    0
sp4    0    0    1    1    1    0
sp5    0    0    0    0    1    1
sp6    0    0    0    0    1    1
> #Numero de especies por local
> S.obs <- apply(cc2,2,sum)
> S.obs
[1] 2 2 3 2 3 2
> ##Maximo de especies em coexistencia
> max(S.obs)
[1] 3
</code>

Podemos nos perguntar se há um limite ao número de espécies que podem coexistir. Se isso ocorre, o número máximo de espécies em co-ocorrência deve ser menor do que se distribuirmos as espécies ao acaso pelos locais. Há vários cenários nulos possíveis, mas ficaremos aqui com permutação dos locais em cada espécie ocorreu. 

Para isso, basta permutar os valores dentro de cada linha:

<code rsplus>
> apply(cc2,1,sample)
     sp1 sp2 sp3 sp4 sp5 sp6
[1,]   1   0   0   1   0   1
[2,]   0   1   0   0   1   0
[3,]   0   1   1   1   0   0
[4,]   0   0   0   1   1   0
[5,]   0   1   0   0   0   1
[6,]   1   0   1   0   0   0
</code>

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 rsplus>
> t(apply(cc2,1,sample))
    [,1] [,2] [,3] [,4] [,5] [,6]
sp1    0    1    0    0    1    0
sp2    1    1    0    0    1    0
sp3    0    1    0    1    0    0
sp4    0    1    0    0    1    1
sp5    1    0    0    0    0    1
sp6    1    0    1    0    0    0
</code>

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 rsplus>
> ##Uma funcao para calcular o maximo de especies em co-ocorrencia
> max.c <- function(x){ max(apply(x,2,sum)) }
> max.coc <- max.c(cc2)
> max.coc
[1] 3
</code>

E então repetimos as permutações 1000 vezes, inserindo-as em um //loop//:

<code rsplus>
##Um vetor para guardar os resultados
> results <- c()
## Mil simulações em um loop
> for(i in 1:1000){
+   result.c[i] <- max.c( t(apply(cc2,1,sample)))
+ }
</code>

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 rsplus>
> sum(result.c<=max.coc)
[1] 216
</code>


===== Reamostragem com Reposição: Bootstrap =====

Como a maioria das boas idéias, o conceito geral do bootstrap é muito simples: 
<box center 80% orange>
//Se aceitamos que nossa amostra é a melhor informação disponível que temos sobre uma população estatística, podemos simular uma nova coleta de dados reamostrando os elementos de nossa amostra, com reposição.//
</box>


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 rsplus>
##Esse objeto de dados está no pacote "vegan"
> data(BCI, package="vegan")
##Visão de três linhas e três colunas
> BCI[1:3,25:28]
  Brosimum.guianense Calophyllum.longifolium Casearia.aculeata Casearia.arborea
1                  0                       0                 0                1
2                  0                       2                 0                1
3                  0                       0                 0                3
</code>

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 rsplus>
> H <- function(x){
+   y <- x[x>0]
+   pi <- y/sum(y)
+   -sum(pi*log(pi))
+ }
> H.BCI <- apply(BCI,1,H)
</code>

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 rsplus>
> hist(H.BCI,xlab="Índice de Shannon")
> qqnorm(H.BCI);qqline(H.BCI)
</code>

<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:
<code rsplus>
> ##Uma amostra ao acaso de 20 parcelas
> bci.sample <- BCI[sample(1:nrow(BCI),20),]

> ##Estimativas da média e intervalo de confiança baseados na distribuição de t de Student:
> t.test(H.am)

        One Sample t-test

data:  H.am
t = 114.8506, df = 19, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 3.783456 3.923914
sample estimates:
mean of x
 3.853685
</code>

A função ''t.test'' quando aplicada a um vetor de números sem nenhum outro argumento retorna o teste t com a hipótese nula de que o vetor é uma amostra de uma população com distribuição normal e média zero. O intervalo de confiança da média da população de onde veio a amostra também é calculado. 

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 rsplus>
> mean(H.BCI)
[1] 3.82084
</code>

Se a amostra vem de uma população de valores com distribuição normal, o intervalo deve incluir a média da população em 95% das vezes que uma amostragem for realizada ((Obs: é um erro comum pensar que o intervalo de confiança é fixo, e a média tem 95% de chance de estar dentro dele. Na verdade, a média populacional é fixa, pois é um parâmetro. O intervalo varia a cada vez que você repetir a amostragem, sendo uma //estimativa por intervalo//)). 

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 rsplus>
##Um vetor para armazenar os resultados
> results <- c()
##Reamostrando 1000 vezes com reposição
> for(i in 1:1000){ results[i] <- mean( sample(H.am,replace=T) )}
</code>

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 rsplus>
> mean(results)
[1] 3.853739
> mean(H.am) - mean(results)
[1] -5.396726e-05
</code>

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 rsplus>
> quantile(results, probs=c(0.025,0.975))
    2.5%    97.5%
3.788261 3.912513
</code>

No caso, não houve diferença importante entre os intervalos de confiança estimados por //bootstrap// e pela distribuição de t de Student.

===== Simulações com Distribuições Teóricas =====
O R gera valores amostrados de distribuições teóricas de probabilidades (ver [[.:03-funcoes#Distribuições Estatísticas: Funções no R]]), o que pode ser usado em simulações de diversas maneiras. 

Um delas é avaliar a eficácia de um procedimento, como a força de um teste de significância em diferentes situações. O exemplo a seguir avalia a força de um teste t comparando duas amostras tomadas de duas populações com distribuição binomial negativa. 

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 rsplus>
> x <- 0:25
> y1 <- dnbinom(x,mu=2, size=2)
> y2 <- dnbinom(x,mu=3, size=2)
> plot(y1~x, type="b", xlab="Valores",ylab="probabilidade",
+      ylim=c(0,max(c(y1,y2))), pch=17 )
> points(y2~x, col="blue", type="b")
> legend(8,0.2,c("média=2, var=4","média=3, var=7,5"),lty=1,
+        col=c("black","blue"), pch=c(17,1))
</code> 

<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**.

Em resumo, usar o teste t para comparar as médias de amostras tomadas de populações com distribuição binomial negativa implica em desconsiderar as premissas do teste, o que pode comprometer sua força. Mas como avaliar isso? Podemos simular amostras tomadas de populações com distribuição binomial negativa com médias diferentes, aplicar o teste t e contar o número de simulações em que o teste indicou uma diferença significativa entre as médias, a 5%.

Para simular a comparação de duas amostras de tamanho dez, usando um //loop// temos:

<code rsplus>
> ##Um vetor fator com 20 elementos, para representar cada elemento das duas amostras
> tratamento <- factor(x=rep(c("a","b"),each=10))
> ##Um vetor para armazenar os resultados da simulação
> results <- c()
> ##Dez mil simulações, com um loop e a função t.test para realizar o teste t
> for(i in 1:10000){
+   variaveis <- c(rnbinom(n=10,size=2,mu=2),rnbinom(n=10,size=2,mu=3))
+   results[i] <- t.test(variaveis~tratamento)$p.value
+ }
> ##Quantas das simulações indicaram uma diferença significativa das médias?
> sum(results<0.05)
[1] 1148
</code>

O teste detectou corretamente que as amostras vêm de populações com médias diferentes em 1148 das dez mil simulações, o que estima a força do teste neste caso em β = 0,1148 ((Este valor irá variar a cada repetição da simulação, pois novos valores serão sorteados, mas ficará em torno do obtido. Para repetir uma simulação e obter os mesmos valores defina a mesma semente de números aleatórios com a função ''set.seed''.)). 

Vamos agora avaliar o que acontece se tomamos amostras de populações com distribuição normal, mantendo as mesmas variâncias:
<code rsplus>
> for(i in 1:10000){
+    variaveis <- c(rnorm(n=10,mean=2,sd=2),
+                   rnorm(n=10,mean=3,sd=sqrt(7.5)))
+    results[i] <- t.test(variaveis~tratamento)$p.value
+  }
> sum(results<0.05)
[1] 1389
</code>

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 rsplus>
> ## Desvio-padrão médio correspondente às variâncias de 7,5 e 4
> sd1 <- mean(sqrt(c(7.5,4)))
> ## Verificando o valor
> sd1
[1] 2.369306
> for(i in 1:10000){
+    variaveis <- c(rnorm(n=10,mean=2,sd=sd1),
+                   rnorm(n=10,mean=3,sd=sd1))
+    results[i] <- t.test(variaveis~tratamento)$p.value
+  }
> sum(results<0.05)
[1] 1419
</code>
 
Mesmo nesta situação ideal com todas as premissas cumpridas, em apenas 14% das tentativas o teste detectará a diferença que existe entre as populações. Isto acontece pelo pequeno tamanho amostral. As simulações podem prosseguir para avaliar o efeito de diferentes tamanhos amostrais.


==== Exercícios ====

<box left 100% | //**Exercício:** Que simulação é esta?//>
Descubra que simulação é feita com código abaixo:
<code rsplus>

> m1 <- matrix(data=rnbinom(n=10000,size=2,mu=2),nrow=10,ncol=1000)
> m2 <- matrix(data=rnbinom(n=10000,size=2,mu=3),nrow=10,ncol=1000)
> matrizona <- rbind(m1,m2)
> tratamento <- factor(x=rep(c("a","b"),each=10))
> signif.t <- function(valores,fator){(t.test(formula=valores~fator))$p.value}
> sum(apply(X=matrizona,MARGIN=2,FUN=signif.t,fator=tratamento)<=0.05)

</code>
</box>