RRR: Reproduzindo o Resampling com Rsampling

Hipótese do estudo

A hipótese deste estudo é que as formigas removem lianas das embaúbas onde estão suas colônias. A previsão é que embaúbas vermelhas seriam menos infestadas por lianas do que as brancas, por abrigarem colônias de formigas mais frequentemente. De fato, esta diferença é observada nas proporções de árvores infestadas na amostra:

> tapply(embauba$with.vines, embauba$morphotype, mean)
   BRANCA  VERMELHA 
0.4878049 0.1000000 

Hipótese nula

A hipótese nula é de que as proporções de infestação são iguais na população de onde vieram as amostras. Sob esta hipótese, uma liana tem a mesma chance de estar em uma embaúba branca ou vermelha. Simulamos a hipótese nula embaralhando as presenças de lianas entre plantas na tabela de dados.

Estatística de interesse

A cada simulação temos que calcular nossa estatística de interesse, que é a a diferença de infestação entre os dois morfos. Criamos uma função para isso:

> emb.ei <- function(dataframe){
+     props <- tapply(dataframe$with.vines, dataframe$morphotype, mean)
+     props[[1]] - props[[2]]
+ }
> ## Verificando
> emb.ei(embauba)
[1] 0.3878049

Distribuição da estatística sob a hipótese nula

Em seguida fazemos a simulação com a função Rsampling:

> emb.r <- Rsampling(type = "normal", dataframe = embauba,
+                    statistics = emb.ei, cols = 2, ntrials = 1000)

O que significa este comando?

• type = "normal" escolhe uma randomização de todos os elementos (mais abaixo você verá outros tipos de randomização).
• dataframe = embauba indica a tabela com os dados
• statistics = emb.ei indica a função que calcula a(s) estatística(s) de interesse da tabela de dados.
• cols = 2 indica que a randomização deve ser feita sobre a segunda coluna da tabela de dados.
• ntrials = 1000 indica o número de repetições da simulação.

A distribuição das estatística de interesse na simulação nem incluiu o valor observado:

> dplot(emb.r, svalue = emb.ei(embauba), pside="Greater",
+       main = "Distribuição da estatística de interesse sob H0",
+       xlab = "Estatística de interesse")

Distribuição das diferenças nas proporções de embaúbas brancas e vermelhas com lianas em 1000 simulações da hipótese nula de ausência de diferença nas populações amostradas. A linha vermelha indica a diferença observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza.

Decisão: rejeitamos a hipótese nula?

Seguindo o padrão das ciências biológicas, adotamos o critério de rejeitar a hipótese nula se a probabilidade da estatística de interesse sob a hipótese nula for menor que 5%.

No gráfico o que não está marcado em cinza são os 5% mais extremos da distribuição da estatística sob a hipótese nula. Então se a estatística observada estiver na região cinza não rejeitamos a hipótese nula. Esta é a chamada de H0. Como o valor observado (linha vermelha) está fora da região de aceitação, podemos rejeitar H0. Você também pode verificar isso com um teste lógico no R:

> sum(emb.r >= emb.ei(embauba))/1000 < 0.05
[1] TRUE

Conclusão: rejeita-se a hipótese nula (p < 0,05).

Hipótese do estudo

A hipótese do estudo é que o recrutamento é mais intenso quando uma folha nova é danificada. A previsão para o experimento é que o recrutamento por extrato de folhas novas seja maior, o que ocorreu:

> splot(azteca$extract.new, azteca$extract.old,
+            group.name=c("Folha nova","Folha velha"),
+            ylab="N de formigas recrutadas",
+            xlab="Tipo de extrato aplicado")

Número de formigas recrutadas por extratos de folhas novas e velhas de embaúbas. Os extratos foram aplicados em pares de folhas próximas em embaúbas que tinham colônias de formigas. As linhas ligam folhas do mesmo par experimental.

Hipótese nula

A hipótese nula é de que o recrutamento provocado pelos estratos é o mesmo. Note que para controlar outras fontes de variação o experimento foi pareado. Então para simular a hipótese nula temos que embaralhar o número de formigas recrutadas dentro de cada par de folhas.

Estatística de interesse

A cada simulação temos que calcular nossa estatística de interesse, que é a média da diferença das folhas de cada par. Uma função para isso:

> azt.ei <- function(dataframe){
+     diferencas <- with(dataframe, extract.new - extract.old)
+     mean(diferencas)
+ }
> ## Valor observado
> azt.ei(azteca)
[1] 7.333333

No experimento o extrato de folhas novas recrutou em média 7.3 formigas que o extrato de folha velha, em cada par.

Distribuição da estatística sob a hipótese nula

Como os pares são linhas em nosso dataframe, simulamos a hipótese nula embaralhando os valores dentro de cada linha:

> azt.r <- Rsampling(type = "within_rows", dataframe = azteca,
+                    statistics = azt.ei, cols = 2:3, ntrials = 1000)

Mudamos o argumento type = "within_rows", para indicar que os valores devem ser embaralhados dentro das linhas. O argumento cols = 2:3 indica as colunas do dataframe que têm as contagens.

Uma diferença igual ou maior que a observada foi muito rara na distribuição da estatística de interesse:

> dplot(azt.r, svalue = azt.ei(azteca), pside="Greater",
+       main = "Distribuição da estatística de interesse sob H0",
+       xlab = "Estatística de interesse")

Distribuição das diferenças do número de formigas recrutadas por extratos de folhas novas e velhas de embaúba em pares experimentais, em 1000 simulações da hipótese nula de ausência de diferença. A linha vermelha indica a diferença observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza.

Decisão: rejeitamos a hipótese nula?

Novamente o gráfico mostra que o valor observado da estatística está fora da região de aceitação da hipótese nula sob nosso critério de significância (5% de chance de erro). O mesmo resultado é verificado com o teste lógico:

> sum(azt.r >= azt.ei(azteca))/1000 < 0.05
[1] TRUE

Conclusão: rejeita-se a hipótese nula (p<0,05).

Coda: testes unicaudais e bicaudais

Até agora testamos hipóteses de que um valor igual ou maior que o observado pode ser gerado pela hipótese nula. É um teste unicaudal ou unidirecional, como seria também o teste de que um valor igual ou menor pode ser gerado pela hipótese nula. Nos testes unicaudais a região de aceitação é toda a distribuição nula exceto seus 5% mais extremos.

Mas pode interessar o teste de que há diferenças, sem especificar sua direção. Por exemplo, o conhecimento prévio poderia indicar a hipótese de que extratos de folhas jovens e velhas devem recrutar números diferentes de formigas, mas sem a expectativa de qual extrato recrutaria mais. Este é um caso de teste bicaudal, quando a região de aceitação é o centro da distribuição nula, exceto seus 2,5% mais extremos de cada lado:

> dplot(azt.r, svalue = azt.ei(azteca), pside="Two sided",
+       main = "Teste bicaudal",
+       xlab = "Estatística de interesse")

Distribuição das diferenças do número de formigas recrutadas por extratos de folhas novas e velhas de embaúba em pares experimentais, em 1000 simulações da hipótese nula de ausência de diferença. A região de aceitação da hipótese nula para 5% de significância para teste bicaudal está delimitada em cinza.

Hipótese do estudo

A hipótese do estudo é que as aranhas preferem caçar em plantas com pelos glandulosos, onde a captura de presas é mais fácil. A previsão para o experimento é que as aranhas devem estar a maior parte do tempo nas folhas com tricomas. De fato, a maioria das aranhas esteve nas folhas com tricomas em 4 ou mais inspeções:

> ## Número de inspeções em que estava em folha com tricomas
> n.insp <- apply(peucetia, 1, sum)
> barplot(table(factor(n.insp, levels=0:6)),
+         xlab="N de inspeções em que estava na folha com tricoma",
+         ylab="N de aranhas")

Número de inspeções em que as 27 aranhas foram registradas em folhas com tricomas, em um experimento de preferência por substratos.

Hipótese nula

A hipótese nula é de que não há preferência. Como metade das placas estavam cobertas com cada tipo de folha, a expectativa nula é que as aranhas estivessem na área coberta por folhas com tricomas em metade das inspeções, em média. Esta expectativa tem a premissa que cada inspeção é um evento independente.

Estatística de interesse

A cada simulação temos que calcular nossa estatística de interesse, que é a média do número de inspeções em que as aranhas estavam sobre folhas com tricomas. Uma função para isso:

> peu.ei <- function(dataframe){
+     mean(apply(dataframe, 1, sum))
+ }
> ## Valor observado
> peu.ei(peucetia)
[1] 4.555556

As aranhas foram registradas em média r round(peu.ei(peucetia),2) das 6 inspeções na área coberta por folhas com tricomas.

Distribuição da estatística sob a hipótese nula

Para simular nossa hipótese nula criamos um data frame com a mesma estrutura, em que cada aranha esteja metade das inspeções nas folhas com tricomas

> peu.H0 <- matrix( rep(c(TRUE,FALSE), each = 3),
+                  nrow = nrow(peucetia), ncol = ncol(peucetia), byrow=TRUE)
> ## Converte em data.frame
> peu.H0 <- data.frame(peu.H0)
> ## verificando
> head(peu.H0)
    X1   X2   X3    X4    X5    X6
1 TRUE TRUE TRUE FALSE FALSE FALSE
2 TRUE TRUE TRUE FALSE FALSE FALSE
3 TRUE TRUE TRUE FALSE FALSE FALSE
4 TRUE TRUE TRUE FALSE FALSE FALSE
5 TRUE TRUE TRUE FALSE FALSE FALSE
6 TRUE TRUE TRUE FALSE FALSE FALSE

E agora simulamos a hipótese nula amostrando com reposição cada linha ²:

> peu.r <- Rsampling(type = "within_rows", dataframe = peu.H0,
+                    statistics = peu.ei, ntrials = 1000, replace=TRUE)

O argumento replace = TRUE, indica amostragem com reposição. No caso isso equivale a sortear uma posição independente para cada aranha a cada inspeção. A probabilidade da aranha estar na folha com tricomas é 0,5 a cada sorteio.

Uma média igual ou maior que a observada não ocorreu na distribuição da estatística de interesse simulada:

> dplot(peu.r, svalue = peu.ei(peucetia), pside="Greater",
+       main = "Distribuição da estatística de interesse sob H0",
+       xlab = "Estatística de interesse")

Distribuição do número médio de inspeções em que as aranhas estavam em folhas com tricomas, em 1000 simulações da hipótese nula de ausência de preferência por substrato. A linha vermelha indica a média observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza.

Decisão: rejeitamos a hipótese nula?

Novamente temos um teste unicaudal, e o valor observado da estatística de interesse não está na região de aceitação da hipótese nula (5%). Confirmamos com o teste lógico do nosso critério de significância:

> sum(peu.r >= peu.ei(peucetia))/1000 < 0.05
[1] TRUE

Conclusão: rejeita-se a hipótese nula (p < 0,05).

Hipótese nula

A proporção das inspeções em que as aranhas permanecem em um dos substratos não depende do tipo de substrato (folha com ou sem tricomas).

Portanto a hipótese nula é sobre a independência entre número de inspeções e tipo de substrato. Simulamos este cenário embaralhando número de ocasiões entre substratos, para cada aranha. Para isso vamos criar um data frame com número de inspeções em cada substrato:

> ## N de inspeções em folha com tricoma
> tric <- apply(peucetia, 1, sum)
> ## N de inspeções em folha lisa
> lisa <- apply(peucetia, 1, function(x) sum(x==0))
> ## Monta o data frame
> peu.H0b <- data.frame(tric=tric, lisa = lisa)
> ## Primeiras linhas
> head(peu.H0b)
  tric lisa
1    0    6
2    0    6
3    0    6
4    2    4
5    3    3
6    3    3

Estatística de interesse

Uma mesma estatística de interesse pode ser aplicada a diferentes hipóteses nulas. Então mantemos a mesma do exemplo anterior: média do número de inspeções em que as aranhas foram registradas nas folhas com tricomas.

Mas como o data frame que será aleatorizado mudou, criamos uma nova função no R para calcular a estatística de interesse

> peu.ei2 <- function(dataframe) mean(dataframe$tric)
> ## Verificando
> peu.ei2(peu.H0b)
[1] 4.555556

Distribuição da estatística de interesse sob a hipótese nula

Simulamos a hipótese nula embaralhando as linhas do data frame com número de inspeções por substrato:

> peu.r2 <- Rsampling(type = "within_rows", dataframe = peu.H0b,
+                    statistics = peu.ei2, ntrials = 1000)

A distribuição nula mudou bastante de forma, comparada com a seção anterior. Mas uma média igual ou maior que a observada ainda foi muito rara:

> dplot(peu.r2, svalue = peu.ei2(peu.H0b), pside="Greater",
+       main = "Distribuição da estatística de interesse sob H0",
+       xlab = "Estatística de interesse")

Distribuição do número médio de inspeções em que as aranhas estavam em folhas com tricomas, em 1000 simulações da hipótese nula de ausência de preferência por substrato, considerando tendência das aranhas permanecerem onde estão. A linha vermelha indica a média observada.

Decisão: rejeitamos a hipótese nula?

O valor observado da estatística de interesse não está na região de aceitação. Aplicando nosso critério de significância:

> sum(peu.r2 >= peu.ei(peucetia))/1000 < 0.05
[1] TRUE

Conclusão: rejeita-se a hipótese nula (p < 0,05).

Hipótese do estudo

Nossa hipótese de estudo é que há ou houve partilha de recursos entre as espécies de pulgões. Neste caso, as associações observadas devem ter resultado em redução da sobreposição de nichos dos insetos.

Hipótese nula

Nossa hipótese nula é que a sobreposição de nicho não difere da esperada caso as ocorrências dos pulgões nas plantas sejam independentes.

Estatística de interesse

Esses dados foram usados para exemplificar um cálculo de amplitude e sobreposição de nichos. A expressão proposta pela autora para a sobreposição média de nichos é a diferença entre a índice de Brillouin de todos os valores e da soma das colunas da tabela. O índice de Brillouin é uma medida de diversidade em uma coleção de valores x_i:

$$H = \frac{1}{N} \log N! \ - \ \frac{1}{N} \sum \log x_i !$$

Onde N = ∑x_i. Vamos criar uma função para fazer este cálculo

> brillouin <- function(x, base=10) {
+     N <- sum(x)
+     lfactorial(N)/(log(base)*N)  -  sum(lfactorial(x)/log(base))/N
+ }

E então criamos um função para o cálculos da estatística de interesse

> pielou.ei <- function(dataframe)
+     brillouin( dataframe ) - brillouin( apply(dataframe,2,sum) )

Cujo valor é

> pielou.ei(pielou)
[1] 0.5804034

Distribuição da estatística de interesse sob a hipótese nula

Para simular nossa hipótese nula, embaralhamos os números de ocorrências registradas de cada espécie de pulgão entre as plantas. Com isso mantemos criamos uma situação em que as hierarquias de preferência de cada espécie de pulgão são mantidas, mas tornam-se independentes. Além disso, usamos a opção fix.zeroes = TRUE para indicar que os valores zero não devem ser embaralhados.

> pielou.r1 <- Rsampling(type = "within_rows", dataframe = pielou,
+                    statistics = pielou.ei, ntrials = 1000, fix.zeroes = TRUE)

O valor observado é maior do que a maioria dos valores na distribuição nula. Como nossa hipótese é unicaudal (sobreposição observada menor do que o esperado pelo acaso), o valor observado está na região de aceitação da hipótese nula.

> dplot(pielou.r1, svalue = pielou.ei(pielou), pside="Lesser",
+       main = "Distribuição da estatística de interesse sob H0",
+       xlab = "Estatística de interesse", xlim=c(0.3,0.6))

Distribuição da sobreposição média de uso de plantas hospedeiras por espécies de pulgões, em 1000 simulações da hipótese nula de independência das espécies de inseto pelas plantas. As plantas sem ocorrência observadas dos pulgões foram consideradas não disponíveis (zeros estruturais). A linha vermelha indica a média observada.

Decisão: rejeitamos a hipótese nula?

O valor observado da estatística de interesse está na região de aceitação. Aplicando nosso critério de significância:

> sum(pielou.r1 <= pielou.ei(pielou))/1000 < 0.05
[1] FALSE

Conclusão: não se rejeita a hipótese nula (p > 0,05).